diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py
index 2ed27c20..e5a63547 100644
--- a/SimPEG/Directives.py
+++ b/SimPEG/Directives.py
@@ -216,7 +216,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the data.
         ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
         phi_ms = 0.5*ms.dot(ms)
-        if self.reg.smoothModel == True:
+        if self.reg.mrefInSmooth == True:
             mref = self.reg.mref
         else:
             mref = 0
@@ -249,7 +249,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the data.
         ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
         phi_ms = 0.5*ms.dot(ms)
-        if self.reg.smoothModel == True:
+        if self.reg.mrefInSmooth == True:
             mref = self.reg.mref
         else:
             mref = 0
@@ -271,7 +271,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
 
-
 # class UpdateReferenceModel(Parameter):
 
 #     mref0 = None
@@ -283,3 +282,63 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 #             mref = self.mref0
 #         self.m_prev = self.invProb.m_current
 #         return mref
+
+class update_IRLS(InversionDirective):
+
+    eps_min = None
+    factor = None
+    gamma = None
+    phi_m_last = None
+    phi_d_last = None
+    
+    def initialize(self):
+
+        # Scale the regularization for changes in norm
+        if getattr(self, 'phi_m_last', None) is not None:
+            self.reg.gamma = 1.
+            phim_new = self.reg.eval(self.invProb.curModel)
+            self.gamma = self.phi_m_last / phim_new
+
+        self.reg.curModel = self.invProb.curModel
+        self.reg.gamma = self.gamma
+        
+        if getattr(self, 'phi_d_last', None) is None:
+            self.phi_d_last = self.invProb.phi_d
+            
+    def endIter(self):
+        # Cool the threshold parameter
+        if getattr(self, 'factor', None) is not None:
+            eps = self.reg.eps / self.factor
+
+            if getattr(self, 'eps_min', None) is not None:
+                self.reg.eps = np.max([self.eps_min,eps])
+            else:
+                self.reg.eps = eps
+
+        # Get phi_m at the end of current iteration
+        self.phi_m_last = self.invProb.phi_m_last
+        
+         # Update the model used for the IRLS weights
+        self.reg.curModel = self.invProb.curModel
+
+        # Update the pre-conditioner
+        diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+        PC     = Utils.sdiag(diagA**-1.)
+        self.opt.approxHinv = PC
+         
+        # Temporarely set gamma to 1.
+        self.reg.gamma = 1.
+        
+        # Compute change in model objective function and update scaling
+        phim_new = self.reg.eval(self.invProb.curModel)
+        self.reg.gamma = self.phi_m_last / phim_new
+        
+        # TO DO: Re-scale beta if too much change in misfit
+        self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d
+        
+#==============================================================================
+#          import pylab as plt
+#          plt.figure()
+#          ax = plt.subplot(221)
+#          self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005))
+#==============================================================================
diff --git a/SimPEG/Examples/MT_1D_ForwardAndInversion.py b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
index cf17e77b..69cedd98 100644
--- a/SimPEG/Examples/MT_1D_ForwardAndInversion.py
+++ b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
@@ -100,7 +100,7 @@ def run(plotIt=True):
     # Regularization - with a regularization mesh
     regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
-    reg.smoothModel = True
+    reg.mrefInSmooth = True
     reg.alpha_s = 1e-7
     reg.alpha_x = 1.
     # Inversion problem
diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py
index 793b2581..d0363001 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')
diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py
index 5b41b91b..c20ed973 100644
--- a/SimPEG/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -1,5 +1,289 @@
 import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
 
+class RegularizationMesh(object):
+    """
+    **Regularization Mesh**
+
+    This contains the operators used in the regularization. Note that these
+    are not necessarily true differential operators, but are constructed from
+    a SimPEG Mesh.
+
+    :param Mesh mesh: problem mesh
+    :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells
+    """
+
+    def __init__(self, mesh, indActive=None):
+        self.mesh = mesh
+        assert indActive is None or indActive.dtype == 'bool', 'indActive needs to be None or a bool'
+        self.indActive = indActive
+
+    @property
+    def vol(self):
+        """
+        reduced volume vector
+        :rtype: numpy.array
+        :return: reduced cell volume
+        """
+        if getattr(self, '_vol', None) is None:
+            self._vol = self._Pac.T * self.mesh.vol
+        return self._vol
+
+    @property
+    def nC(self):
+        """
+        reduced number of cells
+        :rtype: int
+        :return: number of cells being regularized
+        """
+        if getattr(self, '_nC', None) is None:
+            if self.indActive is None:
+                self._nC = self.mesh.nC
+            else:
+                self._nC = sum(self.indActive)
+        return self._nC
+
+    @property
+    def dim(self):
+        """
+        dimension of regularization mesh (1D, 2D, 3D)
+        :rtype: int
+        :return: dimension
+        """
+        if getattr(self, '_dim', None) is None:
+            self._dim = self.mesh.dim
+        return self._dim
+
+
+    @property
+    def _Pac(self):
+        """
+        projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive)
+        :rtype: scipy.sparse.csr_matrix
+        :return: active cell projection matrix
+        """
+        if getattr(self, '__Pac', None) is None:
+            if self.indActive is None:
+                self.__Pac = Utils.speye(self.mesh.nC)
+            else:
+                self.__Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+        return self.__Pac
+
+    @property
+    def _Pafx(self):
+        """
+        projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-x projection matrix
+        """
+        if getattr(self, '__Pafx', None) is None:
+            if self.indActive is None:
+                self.__Pafx = Utils.speye(self.mesh.nFx)
+            else:
+                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
+                self.__Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
+        return self.__Pafx
+
+    @property
+    def _Pafy(self):
+        """
+        projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-y projection matrix
+        """
+        if getattr(self, '__Pafy', None) is None:
+            if self.indActive is None:
+                self.__Pafy = Utils.speye(self.mesh.nFy)
+            else:
+                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
+                self.__Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
+        return self.__Pafy
+
+    @property
+    def _Pafz(self):
+        """
+        projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-z projection matrix
+        """
+        if getattr(self, '__Pafz', None) is None:
+            if self.indActive is None:
+                self.__Pafz = Utils.speye(self.mesh.nFz)
+            else:
+                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):
+        """
+        averaging from active cell centers to active x-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active x-faces
+        """
+        if getattr(self, '_aveFx2CC', None) is None:
+            self._aveFx2CC =  self._Pac.T * self.mesh.aveFx2CC * self._Pafx
+        return self._aveFx2CC
+
+    @property
+    def aveCC2Fx(self):
+        """
+        averaging from active x-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active x-faces to active cell centers
+        """
+        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):
+        """
+        averaging from active cell centers to active y-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active y-faces
+        """
+        if getattr(self, '_aveFy2CC', None) is None:
+            self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy
+        return self._aveFy2CC
+
+    @property
+    def aveCC2Fy(self):
+        """
+        averaging from active y-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active y-faces to active cell centers
+        """
+        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):
+        """
+        averaging from active cell centers to active z-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active z-faces
+        """
+        if getattr(self, '_aveFz2CC', None) is None:
+            self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz
+        return self._aveFz2CC
+
+    @property
+    def aveCC2Fz(self):
+        """
+        averaging from active z-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active z-faces to active cell centers
+        """
+        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 cellDiffx(self):
+        """
+        cell centered difference in the x-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the x-direction
+        """
+        if getattr(self, '_cellDiffx', None) is None:
+            self._cellDiffx = self._Pafx.T * self.mesh.cellGradx * self._Pac
+        return self._cellDiffx
+
+    @property
+    def cellDiffy(self):
+        """
+        cell centered difference in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if getattr(self, '_cellDiffy', None) is None:
+            self._cellDiffy = self._Pafy.T * self.mesh.cellGrady * self._Pac
+        return self._cellDiffy
+
+    @property
+    def cellDiffz(self):
+        """
+        cell centered difference in the z-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the z-direction
+        """
+        if getattr(self, '_cellDiffz', None) is None:
+            self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac
+        return self._cellDiffz
+
+    @property
+    def faceDiffx(self):
+        """
+        x-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the x-direction
+        """
+        if getattr(self, '_faceDiffx', None) is None:
+            self._faceDiffx = self._Pac.T * self.mesh.faceDivx * self._Pafx
+        return self._faceDiffx
+
+    @property
+    def faceDiffy(self):
+        """
+        y-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the y-direction
+        """
+        if getattr(self, '_faceDiffy', None) is None:
+            self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy
+        return self._faceDiffy
+
+    @property
+    def faceDiffz(self):
+        """
+        z-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the z-direction
+        """
+        if getattr(self, '_faceDiffz', None) is None:
+            self._faceDiffz = self._Pac.T * self.mesh.faceDivz * self._Pafz
+        return self._faceDiffz
+
+    @property
+    def cellDiffxStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the x-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the x-direction
+        """
+        if getattr(self, '_cellDiffxStencil', None) is None:
+
+            self._cellDiffxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac
+        return self._cellDiffxStencil
+
+    @property
+    def cellDiffyStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if self.dim < 2: return None
+        if getattr(self, '_cellDiffyStencil', None) is None:
+
+            self._cellDiffyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac
+        return self._cellDiffyStencil
+
+    @property
+    def cellDiffzStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if self.dim < 3: return None
+        if getattr(self, '_cellDiffzStencil', None) is None:
+
+            self._cellDiffzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac
+        return self._cellDiffzStencil
+
+
 class BaseRegularization(object):
     """
     **Base Regularization Class**
@@ -18,12 +302,16 @@ class BaseRegularization(object):
 
     mapping = None    #: A SimPEG.Map instance.
     mesh    = None    #: A SimPEG.Mesh instance.
-    mref = None       #: Reference model.
+    mref    = None    #: Reference model.
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
-        self.mesh = mesh
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        if indActive is not None and indActive.dtype != 'bool':
+            tmp = indActive
+            indActive = np.zeros(mesh.nC, dtype=bool)
+            indActive[tmp] = True
+        self.regmesh = RegularizationMesh(mesh,indActive)
         self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
         self.indActive = indActive
@@ -55,8 +343,8 @@ class BaseRegularization(object):
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.mapping.nP)
-
+        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):
@@ -112,11 +400,10 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
-
 class Tikhonov(BaseRegularization):
     """
     """
-    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    mrefInSmooth = True  #: SMOOTH and SMOOTH_MOD_DIF options
     alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
     alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
     alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
@@ -126,98 +413,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:
-            self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Ws = Pac.T * self._Ws * 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
-            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
-
-            if self.indActive is not None:
-                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
-                self._Wx = Pafx.T*self._Wx*Pac
-
+            Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol
+            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx
         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
-            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
-
-            if self.indActive is not None:
-                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
-                self._Wy = Pafy.T*self._Wy*Pac
-
+            Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol
+            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellDiffy
         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
-            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
-
-            if self.indActive is not None:
-                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
-                self._Wz = Pafz.T*self._Wz*Pac
-
+            Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol
+            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellDiffz
         return self._Wz
 
     @property
     def Wxx(self):
         """Regularization matrix Wxx"""
         if getattr(self, '_Wxx', None) is None:
-            self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wxx = Pac.T*self._Wxx*Pac
-
+            self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDiffx*self.regmesh.cellDiffx
         return self._Wxx
 
     @property
     def Wyy(self):
         """Regularization matrix Wyy"""
         if getattr(self, '_Wyy', None) is None:
-            self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wyy = Pac.T*self._Wyy*Pac
-
+            self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDiffy*self.regmesh.cellDiffy
         return self._Wyy
 
     @property
     def Wzz(self):
         """Regularization matrix Wzz"""
         if getattr(self, '_Wzz', None) is None:
-            self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wzz = Pac.T*self._Wzz*Pac
-
+            self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
         return self._Wzz
 
     @property
@@ -225,9 +472,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
@@ -242,11 +489,11 @@ class Tikhonov(BaseRegularization):
 
     @Utils.timeIt
     def eval(self, m):
-        if self.smoothModel == True:
+        if self.mrefInSmooth == True:
             r1 = self.Wsmooth * ( self.mapping * (m) )
             r2 = self.Ws * ( self.mapping * (m - self.mref) )
             return 0.5*(r1.dot(r1)+r2.dot(r2))
-        elif self.smoothModel == False:
+        elif self.mrefInSmooth == False:
             r = self.W * ( self.mapping * (m - self.mref) )
             return 0.5*r.dot(r)
 
@@ -268,7 +515,7 @@ class Tikhonov(BaseRegularization):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        if self.smoothModel == True:
+        if self.mrefInSmooth == True:
             mD1 = self.mapping.deriv(m)
             mD2 = self.mapping.deriv(m - self.mref)
             r1 = self.Wsmooth * ( self.mapping * (m))
@@ -276,249 +523,223 @@ class Tikhonov(BaseRegularization):
             out1 = mD1.T * ( self.Wsmooth.T * r1 )
             out2 = mD2.T * ( self.Ws.T * r2 )
             out = out1+out2
-        elif self.smoothModel == False:
+        elif self.mrefInSmooth == False:
             mD = self.mapping.deriv(m - self.mref)
             r = self.W * ( self.mapping * (m - self.mref) )
             out = mD.T * ( self.W.T * r )
         return out
 
-# <<<<<<< HEAD
 
-# class Simple(BaseRegularization):
-#     """
-#         Only for tensor mesh
-#     """
+class Simple(BaseRegularization):
+    """
+        Only for tensor mesh
+    """
 
-#     smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
-#     alpha_s  = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
-#     alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-#     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")
-#     alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
-#     alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
-#     alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
+    mrefInSmooth = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
+    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+    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")
+    wght        = 1.
+    
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght
 
-#     def __init__(self, mesh, mapping=None, **kwargs):
-#         BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
+    @property
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self,'_Ws', None) is None:
+                self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
+        return self._Ws
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
+        return self._Wz
+
+    @property
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        if getattr(self, '_Wsmooth', None) is None:
+            wlist = (self.Wx,)
+            if self.regmesh.dim > 1:
+                wlist += (self.Wy,)
+            if self.regmesh.dim > 2:
+                wlist += (self.Wz,)
+            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
+
+
+    @Utils.timeIt
+    def eval(self, m):
+        if self.mrefInSmooth == True:
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
+        elif self.mrefInSmooth == False:
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
+        return phim
 
 
 
-#     @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)
-#         return self._Ws
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
 
-#     @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
-#         return self._Wx
+        The regularization is:
 
-#     @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
-#         return self._Wy
+        .. math::
 
-#     @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
-#         return self._Wz
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
 
-#     @property
-#     def Wxx(self):
-#         """Regularization matrix Wxx"""
-#         if getattr(self, '_Wxx', None) is None:
-#             self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-#         return self._Wxx
+        So the derivative is straight forward:
 
-#     @property
-#     def Wyy(self):
-#         """Regularization matrix Wyy"""
-#         if getattr(self, '_Wyy', None) is None:
-#             self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-#         return self._Wyy
+        .. math::
 
-#     @property
-#     def Wzz(self):
-#         """Regularization matrix Wzz"""
-#         if getattr(self, '_Wzz', None) is None:
-#             self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-#         return self._Wzz
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
-#     @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
-
-#     @Utils.timeIt
-#     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))
-#         elif self.smoothModel == False:
-#             r = self.W * ( self.mapping * (m - self.mref) )
-#             return 0.5*r.dot(r)
+        """
+        if self.mrefInSmooth == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1+out2
+        elif self.mrefInSmooth == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
 
 
-#     @Utils.timeIt
-#     def evalDeriv(self, m):
-#         """
+class Sparse(Simple):
 
-#         The regularization is:
+    # set default values
+    eps      = 1e-1
+    curModel = None # use a model to compute the weights
+    gamma    = 1.
+    p        = 0.
+    qx       = 2.
+    qy       = 2.
+    qz       = 2.
+    wght     = 1.
 
-#         .. math::
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght
 
-#             R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+    @property
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self, 'curModel', None) is None:
+            self.Rs = Utils.speye(self.regmesh.nC)
 
-#         So the derivative is straight forward:
-
-#         .. math::
-
-#             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
-
-#         """
-#         if self.smoothModel == True:
-#             mD1 = self.mapping.deriv(m)
-#             mD2 = self.mapping.deriv(m - self.mref)
-#             r1 = self.Wsmooth * ( self.mapping * (m))
-#             r2 = self.Ws * ( self.mapping * (m - self.mref) )
-#             out1 = mD1.T * ( self.Wsmooth.T * r1 )
-#             out2 = mD2.T * ( self.Ws.T * r2 )
-#             out = out1+out2
-#         elif self.smoothModel == False:
-#             mD = self.mapping.deriv(m - self.mref)
-#             r = self.W * ( self.mapping * (m - self.mref) )
-#             out = mD.T * ( self.W.T * r )
-#         return out
-
-# class SparseRegularization(Simple):
-
-#     eps = 1e-1
-
-#     m = None
-#     gamma = 1.
-#     p = 0.
-#     qx = 2.
-#     qy = 2.
-#     qz = 2.
-
-#     def __init__(self, mesh, mapping=None, **kwargs):
-#         Simple.__init__(self, mesh, mapping=mapping, **kwargs)
+        else:
+            f_m = self.curModel - self.reg.mref
+            self.rs = self.R(f_m , self.p)
+            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
+            self.Rs = Utils.sdiag( self.rs )
+            
+        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
 
 
-#     @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 Wx(self):
+        """Regularization matrix Wx"""
+            
+        if getattr(self, 'curModel', None) is None:
+            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
 
-#     @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
+        else:
+            f_m = self.regmesh.cellDiffxStencil * self.curModel
+            self.rx = self.R( f_m , self.qx)
+            self.Rx = Utils.sdiag( self.rx )
 
-#     @property
-#     def Ws(self):
-#         """Regularization matrix Ws"""
-#         if getattr(self, 'm', None) is None:
-#             self.Rs = Utils.speye(self.mesh.nC)
+        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
 
-#         else:
-#             f_m = self.m
-#             self.rs = self.R(f_m , self.p, self.eps)
-#             #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-#             self.Rs = Utils.sdiag( self.rs )
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
 
-#         self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
+        if getattr(self, 'curModel', None) is None:
+            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
 
-#         return self._Ws
+        else:
+            f_m = self.regmesh.cellDiffyStencil * self.curModel
+            self.ry = self.R( f_m , self.qy)
+            self.Ry = Utils.sdiag( self.ry )
+            
+        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
 
-#     @property
-#     def Wx(self):
-#         """Regularization matrix Wx"""
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
 
-#         if getattr(self, 'm', None) is None:
-#             self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])
+        if getattr(self, 'curModel', None) is None:
+            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
 
-#         else:
-#             f_m = self.mesh.unitCellGradx * self.m
-#             self.rx = self.R( f_m , self.qx, self.eps)
-#             self.Rx = Utils.sdiag( self.rx )
+        else:
+            f_m = self.regmesh.cellDiffzStencil * self.curModel
+            self.rz = self.R( f_m , self.qz)
+            self.Rz = Utils.sdiag( self.rz )
 
-#         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
-#         return self._Wx
+        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
 
-#     @property
-#     def Wy(self):
-#         """Regularization matrix Wy"""
+    @property
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        #if getattr(self, '_Wsmooth', None) is None:
+        wlist = (self.Wx,)
+        if self.regmesh.dim > 1:
+            wlist += (self.Wy,)
+        if self.regmesh.dim > 2:
+            wlist += (self.Wz,)
+        #self._Wsmooth = sp.vstack(wlist)
+        return sp.vstack(wlist)
 
-#         if getattr(self, 'm', None) is None:
-#             self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])
+    @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 sp.vstack(wlist)
+        
+    def R(self, f_m , exponent):
 
-#         else:
-#             f_m = self.mesh.unitCellGrady * self.m
-#             self.ry = self.R( f_m , self.qy, self.eps)
-#             self.Ry = Utils.sdiag( self.ry )
+        eta = (self.eps**(1-exponent/2.))**0.5
+        r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.)
 
-#         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
-#         return self._Wy
-
-#     @property
-#     def Wz(self):
-#         """Regularization matrix Wz"""
-
-#         if getattr(self, 'm', None) is None:
-#             self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])
-
-#         else:
-#             f_m = self.mesh.unitCellGradz * 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
-#         return self._Wz
-
-
-#     def R(self, f_m , p, dec):
-
-#         eta = (self.eps**(1-p/2.))**0.5
-#         r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
-
-#         return r
-# =======
-# >>>>>>> 834de582844e8e1eac95819fbe03eed55dbeb001
+        return r
diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py
index 050c46ac..5aa21ad5 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,44 +18,82 @@ 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):
+
+                    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)
+
+                    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 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
+
+                for i, mesh in enumerate(self.meshlist):
+
+                    print 'Testing Active Cells %iD'%(mesh.dim)
+
+                    if mesh.dim == 1:
+                        indActive = Utils.mkvc(mesh.gridCC <= 0.8)
+                    elif mesh.dim == 2:
+                        indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
+                    elif mesh.dim == 3:
+                        indActive = 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=indActive.nonzero()[0].size)
+
+                    for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
+                        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:', 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)
+
+    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)
-                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:', 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)
-
-    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):
-
-                print 'Testing Active Cells %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)
@@ -62,23 +102,9 @@ class RegularizationTests(unittest.TestCase):
                 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)
+                regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct)
 
-                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:', 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)
+                assert (regmesh.vol == mesh.vol[indAct]).all()
 
 
 if __name__ == '__main__':