diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py
index 85316869..e30bcfaf 100644
--- a/SimPEG/Utils/meshutils.py
+++ b/SimPEG/Utils/meshutils.py
@@ -319,6 +319,83 @@ def writeVTRFile(fileName,mesh,model=None):
     vtrWriteFilter.Update()
 
 
+def ExtractCoremesh(xyzlim, mesh, meshType='tensor'):
+    """
+        Extracts Core Mesh from Global mesh
+        xyzlim: 2D array [ndim x 2]
+        mesh: SimPEG mesh
+        This function ouputs: 
+            - actind: corresponding boolean index from global to core
+            - meshcore: core SimPEG mesh  
+        Warning: 1D and 2D has not been tested
+    """
+    from SimPEG import Mesh
+    if mesh.dim ==1:
+        xyzlim = xyzlim.flatten()
+        xmin, xmax = xyzlim[0], xyzlim[1]
+        
+        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)        
+        
+        xc = mesh.vectorCCx[xind]
+
+        hx = mesh.hx[xind]
+        
+        x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
+        
+        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        
+        actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
+        
+    elif mesh.dim ==2:
+        xmin, xmax = xyzlim[0,0], xyzlim[0,1]
+        ymin, ymax = xyzlim[1,0], xyzlim[1,1]
+
+        yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+
+        xc = mesh.vectorCCx[xind]
+        yc = mesh.vectorCCy[yind]
+
+        hx = mesh.hx[xind]
+        hy = mesh.hy[yind]
+        
+        x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
+        
+        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        
+        actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
+               & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
+        
+    elif mesh.dim==3:
+        xmin, xmax = xyzlim[0,0], xyzlim[0,1]
+        ymin, ymax = xyzlim[1,0], xyzlim[1,1]
+        zmin, zmax = xyzlim[2,0], xyzlim[2,1]
+        
+        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
+        yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+
+        xc = mesh.vectorCCx[xind]
+        yc = mesh.vectorCCy[yind]
+        zc = mesh.vectorCCz[zind]
+
+        hx = mesh.hx[xind]
+        hy = mesh.hy[yind]
+        hz = mesh.hz[zind]
+        
+        x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5, zc[0]-hz[0]*0.5]
+        
+        meshCore = Mesh.TensorMesh([hx, hy, hz] ,x0=x0)
+        
+        actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
+               & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
+               & (mesh.gridCC[:,2]>zmin) & (mesh.gridCC[:,2]<zmax)
+                
+    else:
+        raise(Exception("Not implemented!"))
+    
+    
+    return actind, meshCore
 
 
 if __name__ == '__main__':