diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py
index 22b87edd..10ed2268 100644
--- a/SimPEG/Mesh/DiffOperators.py
+++ b/SimPEG/Mesh/DiffOperators.py
@@ -307,24 +307,28 @@ class DiffOperators(object):
         return BC
     _cellGradBC_list = 'neumann'
 
+    def _cellGradStencil(self):
+        BC = self.setCellGradBC(self._cellGradBC_list)
+        n = self.vnC
+        if(self.dim == 1):
+            G = ddxCellGrad(n[0], BC[0])
+        elif(self.dim == 2):
+            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
+            G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
+            G = sp.vstack((G1, G2), format="csr")
+        elif(self.dim == 3):
+            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
+            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
+            G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            G = sp.vstack((G1, G2, G3), format="csr")
+        return G
+
     def cellGrad():
         doc = "The cell centered Gradient, takes you to cell faces."
 
         def fget(self):
             if(self._cellGrad is None):
-                BC = self.setCellGradBC(self._cellGradBC_list)
-                n = self.vnC
-                if(self.dim == 1):
-                    G = ddxCellGrad(n[0], BC[0])
-                elif(self.dim == 2):
-                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
-                    G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
-                    G = sp.vstack((G1, G2), format="csr")
-                elif(self.dim == 3):
-                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
-                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
-                    G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
-                    G = sp.vstack((G1, G2, G3), format="csr")
+                G = self._cellGradStencil()
                 # Compute areas of cell faces & volumes
                 S = self.area
                 V = self.aveCC2F*self.vol  # Average volume between adjacent cells
@@ -361,19 +365,24 @@ class DiffOperators(object):
     _cellGradBC = None
     cellGradBC = property(**cellGradBC())
 
+    def _cellGradxStencil(self):
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        if(self.dim == 1):
+            G1 = ddxCellGrad(n[0], BC)
+        elif(self.dim == 2):
+            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
+        elif(self.dim == 3):
+            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
+        return G1
+
+
     def cellGradx():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
 
         def fget(self):
             if getattr(self, '_cellGradx', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                if(self.dim == 1):
-                    G1 = ddxCellGrad(n[0], BC)
-                elif(self.dim == 2):
-                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
-                elif(self.dim == 3):
-                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
+                G1 = self._cellGradxStencil()
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fx', 'V')
@@ -382,17 +391,22 @@ class DiffOperators(object):
         return locals()
     cellGradx = property(**cellGradx())
 
+    def _cellGradyStencil(self):
+        if self.dim < 2: return None
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        if(self.dim == 2):
+            G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
+        elif(self.dim == 3):
+            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
+        return G2
+
     def cellGrady():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
             if self.dim < 2: return None
             if getattr(self, '_cellGrady', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                if(self.dim == 2):
-                    G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
-                elif(self.dim == 3):
-                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
+                G2 = self._cellGradyStencil()
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fy', 'V')
@@ -401,14 +415,19 @@ class DiffOperators(object):
         return locals()
     cellGrady = property(**cellGrady())
 
+    def _cellGradzStencil(self):
+        if self.dim < 3: return None
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
+        return G3
+
     def cellGradz():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
             if self.dim < 3: return None
             if getattr(self, '_cellGradz', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
+                G3 = self._cellGradzStencil()
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fz', 'V')
@@ -565,56 +584,6 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
-    def unitCellGradx():
-        doc = """Cell centered Gradient in the x dimension used for
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGradx', None) is None:
-
-                n = self.vnC
-                gx = ddx(n[0]-1)
-                gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr")
-
-                self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square)
-
-            return self._unitCellGradx
-        return locals()
-    unitCellGradx = property(**unitCellGradx())
-
-    def unitCellGrady():
-        doc = """Cell centered Gradient in they dimension used for
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGrady', None) is None:
-
-                n = self.vnC
-                gy = ddx(n[1]-1)
-                gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr")
-
-                self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0]))
-
-            return self._unitCellGrady
-        return locals()
-    unitCellGrady = property(**unitCellGrady())
-
-    def unitCellGradz():
-        doc = """Cell centered Gradient in they dimension used for
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGradz', None) is None:
-
-                n = self.vnC
-                gz = ddx(n[2]-1)
-                gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr")
-
-                self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0]))
-
-            return self._unitCellGradz
-        return locals()
-    unitCellGradz = property(**unitCellGradz())
 
     # --------------- Averaging ---------------------
 
diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py
index d6780c45..31013dc5 100644
--- a/SimPEG/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -1,61 +1,35 @@
 import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
 
-class BaseRegularization(object):
-    """
-    **Base Regularization Class**
+class RegularizationMesh(object):
 
-    This is used to regularize the model space::
-
-        reg = Regularization(mesh)
-
-    """
-
-    __metaclass__ = Utils.SimPEGMetaClass
-
-    counter = None
-
-    mapPair = Maps.IdentityMap    #: A SimPEG.Map Class
-
-    mapping = None    #: A SimPEG.Map instance.
-    mesh    = None    #: A SimPEG.Mesh instance.
-    mref = None       #: Reference model.
-
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        Utils.setKwargs(self, **kwargs)
+    def __init__(self, mesh, indActive=None):
         self.mesh = mesh
-        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
-        self.mapping = mapping or self.mapPair(mesh)
-        self.mapping._assertMatchesPair(self.mapPair)
         self.indActive = indActive
 
     @property
-    def parent(self):
-        """This is the parent of the regularization."""
-        return getattr(self,'_parent',None)
-    @parent.setter
-    def parent(self, p):
-        if getattr(self,'_parent',None) is not None:
-            print 'Regularization has switched to a new parent!'
-        self._parent = p
+    def vol(self):
+        if getattr(self, '_vol', None) is None:
+            self._vol = self._Pac.T * self.mesh.vol
+        return self._vol
 
     @property
-    def inv(self): return self.parent.inv
-    @property
-    def invProb(self): return self.parent
-    @property
-    def reg(self): return self
-    @property
-    def opt(self): return self.parent.opt
-    @property
-    def prob(self): return self.parent.prob
-    @property
-    def survey(self): return self.parent.survey
-
+    def nC(self):
+        if getattr(self, '_nC', None) is None:
+            if self.indActive is None:
+                self._nC = self.mesh.nC
+            else:
+                if self.indActive.dtype == 'bool':
+                    self._nC = sum(self.indActive)
+                else:
+                    self._nC = len(self.indActive) # you shouldn't pass a vector of int 0, 1 's 
+        return self._nC
 
     @property
-    def W(self):
-        """Full regularization weighting matrix W."""
-        return self._Pac.T * sp.identity(self.mesh.nC) * self._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ?
+    def dim(self):
+        if getattr(self, '_dim', None) is None:
+            self._dim = self.mesh.dim
+        return self._dim
+    
 
     @property
     def _Pac(self):
@@ -95,7 +69,170 @@ class BaseRegularization(object):
                 indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
                 self.__Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
         return self.__Pafz
+
+    @property
+    def aveFx2CC(self):
+        if getattr(self, '_aveFx2CC', None) is None:
+            self._aveFx2CC =  self._Pac.T * self.mesh.aveFx2CC * self._Pafx
+        return self._aveFx2CC
+
+    @property
+    def aveCC2Fx(self):
+        if getattr(self, '_aveCC2Fx', None) is None:
+            self._aveCC2Fx =  Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T
+        return self._aveCC2Fx
+
+    @property
+    def aveFy2CC(self):
+        if getattr(self, '_aveFy2CC', None) is None:
+            self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy
+        return self._aveFy2CC
+
+    @property
+    def aveCC2Fy(self):
+        if getattr(self, '_aveCC2Fy', None) is None:
+            self._aveCC2Fy =  Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T
+        return self._aveCC2Fy
+
+    @property
+    def aveFz2CC(self):
+        if getattr(self, '_aveFz2CC', None) is None:
+            self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz
+        return self._aveFz2CC
+
+    @property
+    def aveCC2Fz(self):
+        if getattr(self, '_aveCC2Fz', None) is None:
+            self._aveCC2Fz =  Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T
+        return self._aveCC2Fz
+
+    @property
+    def cellGradx(self):
+        if getattr(self, '_cellGradx', None) is None:
+            self._cellGradx = self._Pafx.T * self.mesh.cellGradx * self._Pac
+        return self._cellGradx
+
+    @property
+    def cellGrady(self):
+        if getattr(self, '_cellGrady', None) is None:
+            self._cellGrady = self._Pafy.T * self.mesh.cellGrady * self._Pac
+        return self._cellGrady
+
+    @property
+    def cellGradz(self):
+        if getattr(self, '_cellGradz', None) is None:
+            self._cellGradz = self._Pafz.T * self.mesh.cellGradz * self._Pac
+        return self._cellGradz
     
+    @property
+    def faceDivx(self):
+        if getattr(self, '_faceDivx', None) is None:
+            self._faceDivx = self._Pac.T * self.mesh.faceDivx * self._Pafx
+        return self._faceDivx
+
+    @property
+    def faceDivy(self):
+        if getattr(self, '_faceDivy', None) is None:
+            self._faceDivy = self._Pac.T * self.mesh.faceDivy * self._Pafy
+        return self._faceDivy
+    
+    @property
+    def faceDivz(self):
+        if getattr(self, '_faceDivz', None) is None:
+            self._faceDivz = self._Pac.T * self.mesh.faceDivz * self._Pafz
+        return self._faceDivz
+
+    @property
+    def cellGradxStencil(self):
+        """Cell centered Gradient in the x dimension used for
+                regularization. The gradient operator is square (nC-by-nC)"""
+
+        # if self.dim < 3: return None
+        if getattr(self, '_cellGradxStencil', None) is None:
+
+            self._cellGradxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac
+        return self._cellGradxStencil
+
+    @property
+    def cellGradyStencil(self):
+        """Cell centered Gradient in the x dimension used for
+                regularization. The gradient operator is square (nC-by-nC)"""
+
+        if self.dim < 2: return None
+        if getattr(self, '_cellGradyStencil', None) is None:
+
+            self._cellGradyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac
+        return self._cellGradyStencil
+
+    @property
+    def cellGradzStencil(self):
+        """Cell centered Gradient in the x dimension used for
+                regularization. The gradient operator is square (nC-by-nC)"""
+
+        if self.dim < 3: return None
+        if getattr(self, '_cellGradzStencil', None) is None:
+
+            self._cellGradzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac
+        return self._cellGradzStencil
+
+
+class BaseRegularization(object):
+    """
+    **Base Regularization Class**
+
+    This is used to regularize the model space::
+
+        reg = Regularization(mesh)
+
+    """
+
+    __metaclass__ = Utils.SimPEGMetaClass
+
+    counter = None
+
+    mapPair = Maps.IdentityMap    #: A SimPEG.Map Class
+
+    mapping = None    #: A SimPEG.Map instance.
+    mesh    = None    #: A SimPEG.Mesh instance.
+    mref    = None    #: Reference model.
+
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        Utils.setKwargs(self, **kwargs)
+        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        self.regmesh = RegularizationMesh(mesh,indActive)
+        self.mapping = mapping or self.mapPair(mesh)
+        self.mapping._assertMatchesPair(self.mapPair)
+        self.indActive = indActive
+
+    @property
+    def parent(self):
+        """This is the parent of the regularization."""
+        return getattr(self,'_parent',None)
+    @parent.setter
+    def parent(self, p):
+        if getattr(self,'_parent',None) is not None:
+            print 'Regularization has switched to a new parent!'
+        self._parent = p
+
+    @property
+    def inv(self): return self.parent.inv
+    @property
+    def invProb(self): return self.parent
+    @property
+    def reg(self): return self
+    @property
+    def opt(self): return self.parent.opt
+    @property
+    def prob(self): return self.parent.prob
+    @property
+    def survey(self): return self.parent.survey
+
+
+    @property
+    def W(self):
+        """Full regularization weighting matrix W."""
+        return sp.identity(self.regmesh.nC)
+        # self.regmesh._Pac.T * sp.identity(self.regmesh.nC) * self.regmesh._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ?
 
     @Utils.timeIt
     def eval(self, m):
@@ -151,7 +288,6 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
-
 class Tikhonov(BaseRegularization):
     """
     """
@@ -165,66 +301,58 @@ class Tikhonov(BaseRegularization):
     alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
 
     def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
-        self.indActive = indActive
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
     @property
     def Ws(self):
         """Regularization matrix Ws"""
         if getattr(self,'_Ws', None) is None:
-            Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-            self._Ws = self._Pac.T * Ws * self._Pac
+            self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
         return self._Ws
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
         if getattr(self, '_Wx', None) is None:
-            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
-            Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx     
-            self._Wx = self._Pafx.T*Wx*self.self._Pac
+            Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol
+            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellGradx     
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
         if getattr(self, '_Wy', None) is None:
-            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
-            Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
-            self._Wy = self._Pafy.T*Wy*self._Pac
+            Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol
+            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellGrady
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
         if getattr(self, '_Wz', None) is None:
-            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
-            Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
-            self._Wz = self._Pafz.T*Wz*self._Pac
+            Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol
+            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellGradz
         return self._Wz
 
     @property
     def Wxx(self):
         """Regularization matrix Wxx"""
         if getattr(self, '_Wxx', None) is None:
-            Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-            self._Wxx = self._Pac.T*Wxx*self._Pac
+            self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDivx*self.regmesh.cellGradx
         return self._Wxx
 
     @property
     def Wyy(self):
         """Regularization matrix Wyy"""
         if getattr(self, '_Wyy', None) is None:
-            Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-            self._Wyy = self._Pac.T*self._Wyy*self._Pac
+            self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDivy*self.regmesh.cellGrady
         return self._Wyy
 
     @property
     def Wzz(self):
         """Regularization matrix Wzz"""
         if getattr(self, '_Wzz', None) is None:
-            Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-            self._Wzz = self._Pac.T*Wzz*self._Pac
+            self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDivz*self.regmesh.cellGradz
         return self._Wzz
 
     @property
@@ -232,9 +360,9 @@ class Tikhonov(BaseRegularization):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
             wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
+            if self.regmesh.dim > 1:
                 wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
+            if self.regmesh.dim > 2:
                 wlist += (self.Wz, self.Wzz)
             self._Wsmooth = sp.vstack(wlist)
         return self._Wsmooth
@@ -301,35 +429,35 @@ class Simple(BaseRegularization):
     alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
     alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
 
-    def __init__(self, mesh, mapping=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
     @property
     def Ws(self):
         """Regularization matrix Ws"""
         if getattr(self,'_Ws', None) is None:
-                self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+                self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
         return self._Ws
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
         if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
+            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x)**0.5)*self.regmesh.cellGradxStencil
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
         if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
+            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y)**0.5)*self.regmesh.cellGradyStencil
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
         if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
+            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z)**0.5)*self.regmesh.cellGradzStencil
         return self._Wz
 
     @property
@@ -337,9 +465,9 @@ class Simple(BaseRegularization):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
             wlist = (self.Wx,)
-            if self.mesh.dim > 1:
+            if self.regmesh.dim > 1:
                 wlist += (self.Wy,)
-            if self.mesh.dim > 2:
+            if self.regmesh.dim > 2:
                 wlist += (self.Wz,)
             self._Wsmooth = sp.vstack(wlist)
         return self._Wsmooth
@@ -352,25 +480,16 @@ class Simple(BaseRegularization):
             self._W = sp.vstack(wlist)
         return self._W
 
-    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.W * ( self.mapping * (m - self.mref) )
-        return 0.5*r.dot(r)
 
     @Utils.timeIt
-    def _evalSmooth(self, m):
+    def eval(self, m):
         if self.smoothModel == True:
             r1 = self.Wsmooth * ( self.mapping * (m) )
             r2 = self.Ws * ( self.mapping * (m - self.mref) )
             return 0.5*(r1.dot(r1)+r2.dot(r2))
-        else:
-            return None
-
-    @Utils.timeIt
-    def eval(self, m):
-        phim = self._evalSmall(m)
-        if self.smoothModel is True:
-            phim += self._evalSmooth(m)
+        elif self.smoothModel == False:
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
         return phim
 
 
@@ -417,34 +536,15 @@ class SparseRegularization(Simple):
     qy    = 2.
     qz    = 2.
 
-    def __init__(self, mesh, mapping=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, **kwargs)
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
 
     @property
     def Ws(self):
         """Regularization matrix Ws"""
         if getattr(self, 'm', None) is None:
-            self.Rs = Utils.speye(self.mesh.nC)
+            self.Rs = Utils.speye(self.regmesh.nC)
 
         else:
             f_m = self.m
@@ -452,7 +552,7 @@ class SparseRegularization(Simple):
             #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
             self.Rs = Utils.sdiag( self.rs )
 
-        self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
+        self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
 
         return self._Ws
 
@@ -461,15 +561,15 @@ class SparseRegularization(Simple):
         """Regularization matrix Wx"""
 
         if getattr(self, 'm', None) is None:
-            self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])
+            self.Rx = Utils.speye(self.regmesh.cellGradxStencil.shape[0])
 
         else:
-            f_m = self.mesh.unitCellGradx * self.m
+            f_m = self.regmesh.cellGradxStencil * self.m
             self.rx = self.R( f_m , self.qx, self.eps)
             self.Rx = Utils.sdiag( self.rx )
 
         if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
+            self._Wx = Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellGradxStencil
         return self._Wx
 
     @property
@@ -477,15 +577,15 @@ class SparseRegularization(Simple):
         """Regularization matrix Wy"""
 
         if getattr(self, 'm', None) is None:
-            self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])
+            self.Ry = Utils.speye(self.regmesh.cellGradyStencil.shape[0])
 
         else:
-            f_m = self.mesh.unitCellGrady * self.m
+            f_m = self.regmesh.cellGradyStencil * self.m
             self.ry = self.R( f_m , self.qy, self.eps)
             self.Ry = Utils.sdiag( self.ry )
 
         if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
+            self._Wy = Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellGradyStencil
         return self._Wy
 
     @property
@@ -493,15 +593,15 @@ class SparseRegularization(Simple):
         """Regularization matrix Wz"""
 
         if getattr(self, 'm', None) is None:
-            self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])
+            self.Rz = Utils.speye(self.regmesh.cellGradzStencil.shape[0])
 
         else:
-            f_m = self.mesh.unitCellGradz * self.m
+            f_m = self.regmesh.cellGradzStencil * self.m
             self.rz = self.R( f_m , self.qz, self.eps)
             self.Rz = Utils.sdiag( self.rz )
 
         if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
+            self._Wz = Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellGradzStencil
         return self._Wz
 
 
diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py
index 050c46ac..614d3158 100644
--- a/tests/base/test_regularization.py
+++ b/tests/base/test_regularization.py
@@ -5,6 +5,8 @@ from scipy.sparse.linalg import dsolve
 import inspect
 
 TOL = 1e-20
+testReg = True
+testRegMesh = True
 
 class RegularizationTests(unittest.TestCase):
 
@@ -16,69 +18,92 @@ class RegularizationTests(unittest.TestCase):
         mesh3 = Mesh.TensorMesh([hx, hy, hz])
         self.meshlist = [mesh1,mesh2, mesh3]
 
-    def test_regularization(self):
-        for R in dir(Regularization):
-            r = getattr(Regularization, R)
-            if not inspect.isclass(r): continue
-            if not issubclass(r, Regularization.BaseRegularization):
-                continue
+    if testReg: 
+        def test_regularization(self):
+            for R in dir(Regularization):
+                r = getattr(Regularization, R)
+                if not inspect.isclass(r): continue
+                if not issubclass(r, Regularization.BaseRegularization):
+                    continue
 
-            for i, mesh in enumerate(self.meshlist):
+                for i, mesh in enumerate(self.meshlist):
 
-                print 'Testing %iD'%mesh.dim
+                    print 'Testing %iD'%mesh.dim
 
-                mapping = r.mapPair(mesh)
-                reg = r(mesh, mapping=mapping)
-                m = np.random.rand(mapping.nP)
-                reg.mref = np.ones_like(m)*np.mean(m)
+                    mapping = r.mapPair(mesh)
+                    reg = r(mesh, mapping=mapping)
+                    m = np.random.rand(mapping.nP)
+                    reg.mref = np.ones_like(m)*np.mean(m)
 
-                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
-                passed = reg.eval(reg.mref) < TOL
-                self.assertTrue(passed)
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                    passed = reg.eval(reg.mref) < TOL
+                    self.assertTrue(passed)
 
-                print 'Check:', R
-                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
+                    print 'Check:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
 
-                print 'Check 2 Deriv:', R
-                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
+                    print 'Check 2 Deriv:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
 
-    def test_regularization_ActiveCells(self):
-        for R in dir(Regularization):
-            r = getattr(Regularization, R)
-            if not inspect.isclass(r): continue
-            if not issubclass(r, Regularization.BaseRegularization):
-                continue
+        def test_regularization_ActiveCells(self):
+            for R in dir(Regularization):
+                r = getattr(Regularization, R)
+                if not inspect.isclass(r): continue
+                if not issubclass(r, Regularization.BaseRegularization):
+                    continue
 
-            for i, mesh in enumerate(self.meshlist):
+                for i, mesh in enumerate(self.meshlist):
 
-                print 'Testing Active Cells %iD'%(mesh.dim)
+                    print 'Testing Active Cells %iD'%(mesh.dim)
 
-                if mesh.dim == 1:
-                    indAct = Utils.mkvc(mesh.gridCC <= 0.8)
-                elif mesh.dim == 2:
-                    indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
-                elif mesh.dim == 3:
-                    indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
+                    if mesh.dim == 1:
+                        indAct = Utils.mkvc(mesh.gridCC <= 0.8)
+                    elif mesh.dim == 2:
+                        indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
+                    elif mesh.dim == 3:
+                        indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
 
-                mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
+                    mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
 
-                reg = r(mesh, mapping=mapping, indActive=indAct)
-                m = np.random.rand(mesh.nC)[indAct]
-                reg.mref = np.ones_like(m)*np.mean(m)
+                    reg = r(mesh, mapping=mapping, indActive=indAct)
+                    m = np.random.rand(mesh.nC)[indAct]
+                    reg.mref = np.ones_like(m)*np.mean(m)
 
-                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
-                passed = reg.eval(reg.mref) < TOL
-                self.assertTrue(passed)
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                    passed = reg.eval(reg.mref) < TOL
+                    self.assertTrue(passed)
 
-                print 'Check:', R
-                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
+                    print 'Check:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
 
-                print 'Check 2 Deriv:', R
-                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
+                    print 'Check 2 Deriv:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+    if testRegMesh:
+        def test_regularizationMesh(self):
+
+                for i, mesh in enumerate(self.meshlist):
+
+                    print 'Testing %iD'%mesh.dim
+
+                    # mapping = r.mapPair(mesh)
+                    # reg = r(mesh, mapping=mapping)
+                    # m = np.random.rand(mapping.nP)
+
+                    if mesh.dim == 1:
+                        indAct = Utils.mkvc(mesh.gridCC <= 0.8)
+                    elif mesh.dim == 2:
+                        indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
+                    elif mesh.dim == 3:
+                        indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
+
+                    regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct)
+
+                    assert (regmesh.vol == mesh.vol[indAct]).all()
 
 
 if __name__ == '__main__':