diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py
index ccde92ae3..1c77f8b88 100644
--- a/discretize/TreeMesh.py
+++ b/discretize/TreeMesh.py
@@ -349,9 +349,10 @@ def permuteE(self):
 
     def plotSlice(
         self, v, vType='CC',
-        normal='Z', ind=None, grid=True, view='real',
+        normal='Z', ind=None, grid=False, view='real',
         ax=None, clim=None, showIt=False,
-        pcolorOpts=None, streamOpts=None, gridOpts=None
+        pcolorOpts=None, streamOpts=None, gridOpts=None,
+        range_x=None, range_y=None,
     ):
 
         if pcolorOpts is None:
@@ -360,10 +361,33 @@ def plotSlice(
             streamOpts = {'color': 'k'}
         if gridOpts is None:
             gridOpts = {'color': 'k', 'alpha': 0.5}
-        if vType not in ['CC', 'F', 'E']:
-            raise ValueError('vType must be one of CC, F, or E')
+        vTypeOpts = ['CC', 'N', 'F', 'E', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', 'Ez']
+        viewOpts = ['real', 'imag', 'abs']
+        normalOpts = ['X', 'Y', 'Z']
+        if vType not in vTypeOpts:
+            raise ValueError(
+                "vType must be in ['{0!s}']".format("', '".join(vTypeOpts))
+            )
         if self.dim == 2:
-            raise Exception('plotSlice not support for 2D, use plotImage')
+            raise NotImplementedError(
+                'Must be a 3D mesh. Use plotImage.'
+            )
+        if view == 'vec':
+            raise NotImplementedError(
+                'Vector view plotting is not implemented for TreeMesh (yet)'
+            )
+        if view not in viewOpts:
+            raise ValueError(
+                "view must be in ['{0!s}']".format("', '".join(viewOpts))
+            )
+        normal = normal.upper()
+        if normal not in normalOpts:
+            raise ValueError(
+                "normal must be in ['{0!s}']".format("', '".join(normalOpts))
+            )
+
+        if not isinstance(grid, bool):
+            raise TypeError('grid must be a boolean')
 
         import matplotlib.pyplot as plt
         import matplotlib
@@ -388,7 +412,7 @@ def plotSlice(
         if type(ind) not in integer_types:
             raise ValueError('ind must be an integer')
 
-        #create a temporary TreeMesh with the slice through
+        # create a temporary TreeMesh with the slice through
         temp_mesh = TreeMesh(h2d, x2d)
         level_diff = self.max_level - temp_mesh.max_level
 
@@ -406,24 +430,41 @@ def plotSlice(
         tm_gridboost[:, antiNormalInd] = temp_mesh.gridCC
         tm_gridboost[:, normalInd] = slice_loc
 
-        #interpolate values to self.gridCC if not 'CC'
-        if vType in ['F', 'E']:
+        # interpolate values to self.gridCC if not 'CC'
+        if vType is not 'CC':
             aveOp = 'ave' + vType + '2CC'
             Av = getattr(self, aveOp)
-            v = Av*v
-
-        #interpolate values from self.gridCC to grid2d
+            if v.size == Av.shape[1]:
+                v = Av*v
+            elif len(vType) == 2:
+                # was one of Fx, Fy, Fz, Ex, Ey, Ez
+                # assuming v has all three components in these cases
+                vec_ind = {'x': 0, 'y': 1, 'z': 2}[vType[1]]
+                if vType[0] == 'E':
+                    i_s = np.cumsum([0, self.nEx, self.nEy, self.nEz])
+                elif vType[0] == 'F':
+                    i_s = np.cumsum([0, self.nFx, self.nFy, self.nFz])
+                v = v[i_s[vec_ind]:i_s[vec_ind+1]]
+                v = Av*v
+
+        # interpolate values from self.gridCC to grid2d
         ind_3d_to_2d = self._get_containing_cell_indexes(tm_gridboost)
         v2d = v[ind_3d_to_2d]
 
         if ax is None:
-            fig = plt.figure()
+            plt.figure()
             ax = plt.subplot(111)
-        else:
-            assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
-            fig = ax.figure
-
-        out = temp_mesh.plotImage(v2d, ax=ax, grid=grid, showIt=showIt, clim=clim)
+        elif not isinstance(ax, matplotlib.axes.Axes):
+            raise Exception("ax must be an matplotlib.axes.Axes")
+
+        out = temp_mesh.plotImage(
+            v2d, vType='CC',
+            grid=grid, view=view,
+            ax=ax, clim=clim, showIt=False,
+            pcolorOpts=pcolorOpts,
+            gridOpts=gridOpts,
+            range_x=range_x,
+            range_y=range_y)
 
         ax.set_xlabel('y' if normal == 'X' else 'x')
         ax.set_ylabel('y' if normal == 'Z' else 'z')
diff --git a/discretize/tree.cpp b/discretize/tree.cpp
index 775c23479..93b7d16d0 100644
--- a/discretize/tree.cpp
+++ b/discretize/tree.cpp
@@ -395,148 +395,160 @@ void Cell::divide(node_map_t& nodes, double* xs, double* ys, double* zs, bool fo
     }else if(force){
         do_splitting = true;
     }else{
-        int_t test_level = (*test_func)(this);
-        if(test_level > level){
-            do_splitting = true;
+        if(is_leaf()){
+            // Only need to call the function if I am a leaf...
+            int_t test_level = (*test_func)(this);
+            do_splitting = test_level > level;
         }
     }
-    if(!do_splitting){
-        return;
-    }
-    //If i haven't already been split...
-    if(children[0] == NULL){
-        spawn(nodes, children, xs, ys, zs);
-
-        //If I need to be split, and my neighbor is below my level
-        //Then it needs to be split
-        //-x,+x,-y,+y,-z,+z
-        if(balance){
-            for(int_t i = 0; i < 2*n_dim; ++i){
-                if(neighbors[i] != NULL && neighbors[i]->level < level){
-                    neighbors[i]->divide(nodes, xs, ys, zs, true);
+    if(do_splitting){
+        //If i haven't already been split...
+        if(is_leaf()){
+            spawn(nodes, children, xs, ys, zs);
+
+            //If I need to be split, and my neighbor is below my level
+            //Then it needs to be split
+            //-x,+x,-y,+y,-z,+z
+            if(balance){
+                for(int_t i = 0; i < 2*n_dim; ++i){
+                    if(neighbors[i] != NULL && neighbors[i]->level < level){
+                        neighbors[i]->divide(nodes, xs, ys, zs, true);
+                    }
                 }
             }
-        }
 
-        //Set children's neighbors (first do the easy ones)
-        // all of the children live next to each other
-        children[0]->set_neighbor(children[1], 1);
-        children[0]->set_neighbor(children[2], 3);
-        children[1]->set_neighbor(children[3], 3);
-        children[2]->set_neighbor(children[3], 1);
-
-        if(n_dim == 3){
-            children[4]->set_neighbor(children[5], 1);
-            children[4]->set_neighbor(children[6], 3);
-            children[5]->set_neighbor(children[7], 3);
-            children[6]->set_neighbor(children[7], 1);
-
-            children[0]->set_neighbor(children[4], 5);
-            children[1]->set_neighbor(children[5], 5);
-            children[2]->set_neighbor(children[6], 5);
-            children[3]->set_neighbor(children[7], 5);
-        }
+            //Set children's neighbors (first do the easy ones)
+            // all of the children live next to each other
+            children[0]->set_neighbor(children[1], 1);
+            children[0]->set_neighbor(children[2], 3);
+            children[1]->set_neighbor(children[3], 3);
+            children[2]->set_neighbor(children[3], 1);
+
+            if(n_dim == 3){
+                children[4]->set_neighbor(children[5], 1);
+                children[4]->set_neighbor(children[6], 3);
+                children[5]->set_neighbor(children[7], 3);
+                children[6]->set_neighbor(children[7], 1);
+
+                children[0]->set_neighbor(children[4], 5);
+                children[1]->set_neighbor(children[5], 5);
+                children[2]->set_neighbor(children[6], 5);
+                children[3]->set_neighbor(children[7], 5);
+            }
 
-        // -x direction
-        if(neighbors[0]!=NULL && !(neighbors[0]->is_leaf())){
-            children[0]->set_neighbor(neighbors[0]->children[1], 0);
-            children[2]->set_neighbor(neighbors[0]->children[3], 0);
-        }
-        else{
-            children[0]->set_neighbor(neighbors[0], 0);
-            children[2]->set_neighbor(neighbors[0], 0);
-        }
-        // +x direction
-        if(neighbors[1]!=NULL && !neighbors[1]->is_leaf()){
-            children[1]->set_neighbor(neighbors[1]->children[0], 1);
-            children[3]->set_neighbor(neighbors[1]->children[2], 1);
-        }else{
-            children[1]->set_neighbor(neighbors[1], 1);
-            children[3]->set_neighbor(neighbors[1], 1);
-        }
-        // -y direction
-        if(neighbors[2]!=NULL && !neighbors[2]->is_leaf()){
-            children[0]->set_neighbor(neighbors[2]->children[2], 2);
-            children[1]->set_neighbor(neighbors[2]->children[3], 2);
-        }else{
-            children[0]->set_neighbor(neighbors[2], 2);
-            children[1]->set_neighbor(neighbors[2], 2);
-        }
-        // +y direction
-        if(neighbors[3]!=NULL && !neighbors[3]->is_leaf()){
-            children[2]->set_neighbor(neighbors[3]->children[0], 3);
-            children[3]->set_neighbor(neighbors[3]->children[1], 3);
-        }else{
-            children[2]->set_neighbor(neighbors[3], 3);
-            children[3]->set_neighbor(neighbors[3], 3);
-        }
-        if(n_dim==3){
             // -x direction
             if(neighbors[0]!=NULL && !(neighbors[0]->is_leaf())){
-                children[4]->set_neighbor(neighbors[0]->children[5], 0);
-                children[6]->set_neighbor(neighbors[0]->children[7], 0);
+                children[0]->set_neighbor(neighbors[0]->children[1], 0);
+                children[2]->set_neighbor(neighbors[0]->children[3], 0);
             }
             else{
-                children[4]->set_neighbor(neighbors[0], 0);
-                children[6]->set_neighbor(neighbors[0], 0);
+                children[0]->set_neighbor(neighbors[0], 0);
+                children[2]->set_neighbor(neighbors[0], 0);
             }
             // +x direction
             if(neighbors[1]!=NULL && !neighbors[1]->is_leaf()){
-                children[5]->set_neighbor(neighbors[1]->children[4], 1);
-                children[7]->set_neighbor(neighbors[1]->children[6], 1);
+                children[1]->set_neighbor(neighbors[1]->children[0], 1);
+                children[3]->set_neighbor(neighbors[1]->children[2], 1);
             }else{
-                children[5]->set_neighbor(neighbors[1], 1);
-                children[7]->set_neighbor(neighbors[1], 1);
+                children[1]->set_neighbor(neighbors[1], 1);
+                children[3]->set_neighbor(neighbors[1], 1);
             }
             // -y direction
             if(neighbors[2]!=NULL && !neighbors[2]->is_leaf()){
-                children[4]->set_neighbor(neighbors[2]->children[6], 2);
-                children[5]->set_neighbor(neighbors[2]->children[7], 2);
+                children[0]->set_neighbor(neighbors[2]->children[2], 2);
+                children[1]->set_neighbor(neighbors[2]->children[3], 2);
             }else{
-                children[4]->set_neighbor(neighbors[2], 2);
-                children[5]->set_neighbor(neighbors[2], 2);
+                children[0]->set_neighbor(neighbors[2], 2);
+                children[1]->set_neighbor(neighbors[2], 2);
             }
             // +y direction
             if(neighbors[3]!=NULL && !neighbors[3]->is_leaf()){
-                children[6]->set_neighbor(neighbors[3]->children[4], 3);
-                children[7]->set_neighbor(neighbors[3]->children[5], 3);
+                children[2]->set_neighbor(neighbors[3]->children[0], 3);
+                children[3]->set_neighbor(neighbors[3]->children[1], 3);
             }else{
-                children[6]->set_neighbor(neighbors[3], 3);
-                children[7]->set_neighbor(neighbors[3], 3);
+                children[2]->set_neighbor(neighbors[3], 3);
+                children[3]->set_neighbor(neighbors[3], 3);
             }
-            // -z direction
-            if(neighbors[4]!=NULL && !neighbors[4]->is_leaf()){
-                children[0]->set_neighbor(neighbors[4]->children[4], 4);
-                children[1]->set_neighbor(neighbors[4]->children[5], 4);
-                children[2]->set_neighbor(neighbors[4]->children[6], 4);
-                children[3]->set_neighbor(neighbors[4]->children[7], 4);
-            }else{
-                children[0]->set_neighbor(neighbors[4], 4);
-                children[1]->set_neighbor(neighbors[4], 4);
-                children[2]->set_neighbor(neighbors[4], 4);
-                children[3]->set_neighbor(neighbors[4], 4);
-            }
-            // +z direction
-            if(neighbors[5]!=NULL && !neighbors[5]->is_leaf()){
-                children[4]->set_neighbor(neighbors[5]->children[0], 5);
-                children[5]->set_neighbor(neighbors[5]->children[1], 5);
-                children[6]->set_neighbor(neighbors[5]->children[2], 5);
-                children[7]->set_neighbor(neighbors[5]->children[3], 5);
-            }else{
-                children[4]->set_neighbor(neighbors[5], 5);
-                children[5]->set_neighbor(neighbors[5], 5);
-                children[6]->set_neighbor(neighbors[5], 5);
-                children[7]->set_neighbor(neighbors[5], 5);
+            if(n_dim==3){
+                // -x direction
+                if(neighbors[0]!=NULL && !(neighbors[0]->is_leaf())){
+                    children[4]->set_neighbor(neighbors[0]->children[5], 0);
+                    children[6]->set_neighbor(neighbors[0]->children[7], 0);
+                }
+                else{
+                    children[4]->set_neighbor(neighbors[0], 0);
+                    children[6]->set_neighbor(neighbors[0], 0);
+                }
+                // +x direction
+                if(neighbors[1]!=NULL && !neighbors[1]->is_leaf()){
+                    children[5]->set_neighbor(neighbors[1]->children[4], 1);
+                    children[7]->set_neighbor(neighbors[1]->children[6], 1);
+                }else{
+                    children[5]->set_neighbor(neighbors[1], 1);
+                    children[7]->set_neighbor(neighbors[1], 1);
+                }
+                // -y direction
+                if(neighbors[2]!=NULL && !neighbors[2]->is_leaf()){
+                    children[4]->set_neighbor(neighbors[2]->children[6], 2);
+                    children[5]->set_neighbor(neighbors[2]->children[7], 2);
+                }else{
+                    children[4]->set_neighbor(neighbors[2], 2);
+                    children[5]->set_neighbor(neighbors[2], 2);
+                }
+                // +y direction
+                if(neighbors[3]!=NULL && !neighbors[3]->is_leaf()){
+                    children[6]->set_neighbor(neighbors[3]->children[4], 3);
+                    children[7]->set_neighbor(neighbors[3]->children[5], 3);
+                }else{
+                    children[6]->set_neighbor(neighbors[3], 3);
+                    children[7]->set_neighbor(neighbors[3], 3);
+                }
+                // -z direction
+                if(neighbors[4]!=NULL && !neighbors[4]->is_leaf()){
+                    children[0]->set_neighbor(neighbors[4]->children[4], 4);
+                    children[1]->set_neighbor(neighbors[4]->children[5], 4);
+                    children[2]->set_neighbor(neighbors[4]->children[6], 4);
+                    children[3]->set_neighbor(neighbors[4]->children[7], 4);
+                }else{
+                    children[0]->set_neighbor(neighbors[4], 4);
+                    children[1]->set_neighbor(neighbors[4], 4);
+                    children[2]->set_neighbor(neighbors[4], 4);
+                    children[3]->set_neighbor(neighbors[4], 4);
+                }
+                // +z direction
+                if(neighbors[5]!=NULL && !neighbors[5]->is_leaf()){
+                    children[4]->set_neighbor(neighbors[5]->children[0], 5);
+                    children[5]->set_neighbor(neighbors[5]->children[1], 5);
+                    children[6]->set_neighbor(neighbors[5]->children[2], 5);
+                    children[7]->set_neighbor(neighbors[5]->children[3], 5);
+                }else{
+                    children[4]->set_neighbor(neighbors[5], 5);
+                    children[5]->set_neighbor(neighbors[5], 5);
+                    children[6]->set_neighbor(neighbors[5], 5);
+                    children[7]->set_neighbor(neighbors[5], 5);
+                }
             }
         }
     }
     if(!force){
-        for(int_t i = 0; i < (1<<n_dim); ++i){
-            children[i]->divide(nodes, xs, ys, zs);
+        if(!is_leaf()){
+            for(int_t i = 0; i < (1<<n_dim); ++i){
+                children[i]->divide(nodes, xs, ys, zs);
+            }
         }
     }
 };
 
+
+void Cell::set_test_function(function func){
+    test_func = func;
+    if(!is_leaf()){
+        for(int_t i = 0; i < (1<<n_dim); ++i){
+            children[i]->set_test_function(func);
+        }
+    }
+}
+
 void Cell::build_cell_vector(cell_vec_t& cells){
     if(this->is_leaf()){
         cells.push_back(this);
@@ -723,7 +735,7 @@ void Tree::build_tree_from_function(function test_func){
     for(int_t iz=0; iz<nz_roots; ++iz)
         for(int_t iy=0; iy<ny_roots; ++iy)
             for(int_t ix=0; ix<nx_roots; ++ix)
-                roots[iz][iy][ix]->test_func = test_func;
+                roots[iz][iy][ix]->set_test_function(test_func);
     //Now we can divide
     for(int_t iz=0; iz<nz_roots; ++iz)
         for(int_t iy=0; iy<ny_roots; ++iy)
diff --git a/discretize/tree.h b/discretize/tree.h
index b4e98b6ae..106977855 100644
--- a/discretize/tree.h
+++ b/discretize/tree.h
@@ -120,6 +120,7 @@ class Cell{
     void spawn(node_map_t& nodes, Cell *kids[8], double* xs, double *ys, double *zs);
     void divide(node_map_t& nodes, double* xs, double* ys, double* zs, bool force=false, bool balance=true);
     void set_neighbor(Cell* other, int_t direction);
+    void set_test_function(function func);
     void build_cell_vector(cell_vec_t& cells);
 
     void insert_cell(node_map_t &nodes, double *new_center, int_t p_level, double* xs, double *ys, double *zs);
diff --git a/discretize/tree_ext.pyx b/discretize/tree_ext.pyx
index c5b974fb2..74c216ca7 100644
--- a/discretize/tree_ext.pyx
+++ b/discretize/tree_ext.pyx
@@ -144,6 +144,7 @@ cdef class _TreeMesh:
     cdef PyWrapper *wrapper
     cdef int_t _dim
     cdef int_t[3] ls
+    cdef int _finalized
 
     cdef double[:] _xs, _ys, _zs
     cdef double[:] _x0
@@ -202,7 +203,10 @@ cdef class _TreeMesh:
         self.tree.set_levels(self.ls[0], self.ls[1], self.ls[2])
         self.tree.set_xs(&self._xs[0], &self._ys[0], &self._zs[0])
         self.tree.initialize_roots()
+        self._finalized = False
+        self._clear_cache()
 
+    def _clear_cache(self):
         self._gridCC = None
         self._gridN = None
         self._gridhN = None
@@ -277,7 +281,8 @@ cdef class _TreeMesh:
         points = np.require(np.atleast_2d(points), dtype=np.float64,
                             requirements='C')
         cdef double[:, :] cs = points
-        cdef int[:] ls = np.asarray(levels, dtype=np.int32)
+        cdef int[:] ls = np.require(np.atleast_1d(levels), dtype=np.int32,
+                                    requirements='C')
         cdef int_t i
         for i in range(ls.shape[0]):
             self.tree.insert_cell(&cs[i, 0], ls[i])
@@ -285,8 +290,10 @@ cdef class _TreeMesh:
             self.finalize()
 
     def finalize(self):
-        self.tree.finalize_lists()
-        self.tree.number()
+        if not self._finalized:
+            self.tree.finalize_lists()
+            self.tree.number()
+            self._finalized=True
 
     def number(self):
         self.tree.number()
@@ -3012,9 +3019,22 @@ cdef class _TreeMesh:
         if showIt:
             plt.show()
 
-    def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
+    def plotImage(self, v, vType='CC', grid=False, view='real',
+                  ax=None, clim=None, showIt=False,
+                  pcolorOpts=None,
+                  gridOpts=None,
+                  range_x=None, range_y=None,
+                  **other_kwargs,
+                  ):
         if self._dim == 3:
-            raise Exception('Use plot slice')
+            self.plotSlice(v, vType=vType, grid=grid, view=view,
+                           ax=ax, clim=clim, showIt=showIt,
+                           pcolorOpts=pcolorOpts,
+                           range_x=range_x, range_y=range_y,
+                           **other_kwargs)
+
+        if view == 'vec':
+            raise NotImplementedError('Vector ploting is not supported on TreeMesh (yet)')
 
         import matplotlib.pyplot as plt
         import matplotlib
@@ -3022,19 +3042,65 @@ cdef class _TreeMesh:
         import matplotlib.cm as cmx
         from matplotlib.collections import PatchCollection
 
+        if vType == 'CC':
+            I = v
+        elif vType == 'N':
+            I = self.aveN2CC*v
+        elif vType in ['Fx', 'Fy', 'Ex', 'Ey']:
+            aveOp = 'ave' + vType[0] + '2CCV'
+            ind_xy = {'x': 0, 'y': 1}[vType[1]]
+            I = (getattr(self, aveOp)*v).reshape(2, self.nC)[ind_xy] # average to cell centers
+
+        if view in ['real', 'imag', 'abs']:
+            v = getattr(np, view)(v) # e.g. np.real(v)
         if ax is None:
             ax = plt.subplot(111)
-        default = cm = plt.get_cmap()
+        if pcolorOpts is None:
+            pcolorOpts = {}
+        if 'cmap' in pcolorOpts:
+            cm = pcolorOpts['cmap']
+        else:
+            cm = plt.get_cmap()
+        if 'vmin' in pcolorOpts:
+            vmin = pcolorOpts['vmin']
+        else:
+            vmin = np.nanmin(I) if clim is None else clim[0]
+        if 'vmax' in pcolorOpts:
+            vmax = pcolorOpts['vmax']
+        else:
+            vmax = np.nanmax(I) if clim is None else clim[1]
+        if 'alpha' in pcolorOpts:
+            alpha = pcolorOpts['alpha']
+        else:
+            alpha = 1.0
+
+        if gridOpts is None:
+            gridOpts = {'color':'k'}
+        if 'color' not in gridOpts:
+            gridOpts['color'] = 'k'
+
         cNorm = colors.Normalize(
-            vmin=np.nanmin(I) if clim is None else clim[0],
-            vmax=np.nanmax(I) if clim is None else clim[1])
+            vmin=vmin, vmax=vmax)
+
+        scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
+
+        if 'edge_color' in pcolorOpts:
+            edge_color = pcolorOpts['edge_color']
+        else:
+            edge_color = gridOpts['color'] if grid else 'none'
+        if 'alpha' in gridOpts:
+            edge_alpha = gridOpts['alpha']
+        else:
+            edge_alpha = 1.0
+        if edge_color.lower() != 'none':
+            if isinstance(edge_color, str):
+                edge_color = colors.to_rgba(edge_color, edge_alpha)
+            else:
+                edge_color = colors.to_rgba_array(edge_color, edge_alpha)
 
-        scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=default)
-        ax.set_xlim((self.x0[0], self.h[0].sum()))
-        ax.set_ylim((self.x0[1], self.h[1].sum()))
-        edge_color = 'k' if grid else 'none'
         rectangles = []
-        facecolors = []
+        facecolors = scalarMap.to_rgba(I[~np.isnan(I)])
+        facecolors[:,-1] = alpha
         for cell in self.tree.cells:
             ii = cell.index
             if np.isnan(I[ii]):
@@ -3042,7 +3108,6 @@ cdef class _TreeMesh:
             x0 = np.array([cell.points[0].location[0], cell.points[0].location[1]])
             sz = np.array([cell.edges[0].length, cell.edges[2].length])
             rectangles.append(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1]))
-            facecolors.append(scalarMap.to_rgba(I[ii]))
 
         pc = PatchCollection(rectangles, facecolor=facecolors, edgecolor=edge_color)
         # Add collection to axes
@@ -3051,6 +3116,17 @@ cdef class _TreeMesh:
         scalarMap._A = []
         ax.set_xlabel('x')
         ax.set_ylabel('y')
+
+        if range_x is not None:
+            ax.set_xlim(*range_x)
+        else:
+            ax.set_xlim(*self.vectorNx[[0, -1]])
+
+        if range_y is not None:
+            ax.set_ylim(*range_y)
+        else:
+            ax.set_ylim(*self.vectorNy[[0, -1]])
+
         if showIt:
             plt.show()
         return [scalarMap]
diff --git a/discretize/utils/__init__.py b/discretize/utils/__init__.py
index 34d78e983..22b4d7109 100644
--- a/discretize/utils/__init__.py
+++ b/discretize/utils/__init__.py
@@ -10,7 +10,7 @@
 from .codeutils import (isScalar, asArray_N_x_Dim)
 from .meshutils import (
     exampleLrmGrid, meshTensor, closestPoints, ExtractCoreMesh,
-    random_model
+    random_model, mesh_builder_xyz, refine_tree_xyz
 )
 from .curvutils import volTetra, faceInfo, indexCube
 from .interputils import interpmat
diff --git a/discretize/utils/meshutils.py b/discretize/utils/meshutils.py
index 19c1ddcf9..1c3830ff4 100644
--- a/discretize/utils/meshutils.py
+++ b/discretize/utils/meshutils.py
@@ -5,6 +5,9 @@
 from .matutils import ndgrid
 from .codeutils import asArray_N_x_Dim
 from .codeutils import isScalar
+import discretize
+from scipy.spatial import cKDTree, Delaunay
+from scipy import interpolate
 
 import sys
 if sys.version_info < (3,):
@@ -38,25 +41,38 @@ def exampleLrmGrid(nC, exType):
 
 def random_model(shape, seed=None, anisotropy=None, its=100, bounds=None):
     """
-        Create a random model by convolving a kernel with a
-        uniformly distributed model.
+    Create a random model by convolving a kernel with a
+    uniformly distributed model.
+
+    Parameters
+    ----------
+    shape: tuple
+        shape of the model.
+    seed: int
+        pick which model to produce, prints the seed if you don't choose.
+    anisotropy: numpy.ndarray
+        this is the (3 x n) blurring kernel that is used.
+    its: int
+        number of smoothing iterations
+    bounds: list
+        bounds on the model, len(list) == 2
+
+    Returns
+    -------
+    numpy.ndarray
+        M, the model
+
+    Example
+    -------
 
-        :param tuple shape: shape of the model.
-        :param int seed: pick which model to produce, prints the seed if you don't choose.
-        :param numpy.ndarray anisotropy: this is the (3 x n) blurring kernel that is used.
-        :param int its: number of smoothing iterations
-        :param list bounds: bounds on the model, len(list) == 2
-        :rtype: numpy.ndarray
-        :return: M, the model
-
-
-        .. plot::
+    .. plot::
+        :include-source:
 
-            import matplotlib.pyplot as plt
-            import discretize
-            plt.colorbar(plt.imshow(discretize.utils.random_model((50, 50), bounds=[-4, 0])))
-            plt.title('A very cool, yet completely random model.')
-            plt.show()
+        import matplotlib.pyplot as plt
+        import discretize
+        plt.colorbar(plt.imshow(discretize.utils.random_model((50, 50), bounds=[-4, 0])))
+        plt.title('A very cool, yet completely random model.')
+        plt.show()
 
 
     """
@@ -96,7 +112,6 @@ def random_model(shape, seed=None, anisotropy=None, its=100, bounds=None):
     return mi
 
 
-
 def meshTensor(value):
     """**meshTensor** takes a list of numbers and tuples
     that have the form::
@@ -120,12 +135,13 @@ def meshTensor(value):
     is negative this section of the tensor is flipped right-to-left.
 
     .. plot::
+        :include-source:
 
         import discretize
         tx = [(10.0, 10, -1.3), (10.0, 40), (10.0, 10, 1.3)]
         ty = [(10.0, 10, -1.3), (10.0, 40)]
-        M = discretize.TensorMesh([tx, ty])
-        M.plotGrid(showIt=True)
+        mesh = discretize.TensorMesh([tx, ty])
+        mesh.plotGrid(showIt=True)
 
     """
     if type(value) is not list:
@@ -153,11 +169,20 @@ def meshTensor(value):
 def closestPoints(mesh, pts, gridLoc='CC'):
     """Move a list of points to the closest points on a grid.
 
-    :param BaseMesh mesh: The mesh
-    :param numpy.ndarray pts: Points to move
-    :param str gridLoc: ['CC', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ex', 'Ey', 'Ez']
-    :rtype: numpy.ndarray
-    :return: nodeInds
+    Parameters
+    ----------
+    mesh: BaseMesh
+        The mesh
+    pts: numpy.ndarray
+        Points to move
+    gridLoc: str
+        ['CC', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ex', 'Ey', 'Ez']
+
+    Returns
+    -------
+    numpy.ndarray
+        nodeInds
+
     """
 
     pts = asArray_N_x_Dim(pts, mesh.dim)
@@ -173,19 +198,25 @@ def closestPoints(mesh, pts, gridLoc='CC'):
     return nodeInds
 
 
-def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
+def ExtractCoreMesh(xyzlim, mesh, mesh_type='tensor'):
     """
     Extracts Core Mesh from Global mesh
 
-    :param numpy.ndarray xyzlim: 2D array [ndim x 2]
-    :param BaseMesh mesh: The mesh
-
-    This function ouputs::
-
-        - actind: corresponding boolean index from global to core
-        - meshcore: core sdiscretize mesh
+    Parameters
+    ----------
+    xyzlim: numpy.ndarray
+        2D array [ndim x 2]
+    mesh: BaseMesh
+        The mesh
+
+    Returns
+    -------
+    tuple
+        This function ouputs::
+            - actind: corresponding boolean index from global to core
+            - meshcore: core sdiscretize mesh
     """
-    import discretize
+
     if mesh.dim == 1:
         xyzlim = xyzlim.flatten()
         xmin, xmax = xyzlim[0], xyzlim[1]
@@ -255,3 +286,487 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
         raise Exception("Not implemented!")
 
     return actind, meshCore
+
+
+def mesh_builder_xyz(
+    xyz, h,
+    padding_distance=[[0, 0], [0, 0], [0, 0]],
+    base_mesh=None,
+    depth_core=None,
+    expansion_factor=1.3,
+    mesh_type='tensor'
+):
+    """
+    Function to quickly generate a Tensor or Tree mesh
+    given a cloud of xyz points, finest core cell size
+    and padding distance.
+    If a base_mesh is provided, the core cells will be centered
+    on the underlaying mesh to reduce interpolation errors.
+    The core extent is set by the bounding box of the xyz location.
+
+    Parameters
+    ----------
+    xyz: numpy.ndarray
+        Location points [n x dim]
+    h: list
+        Cell size for the core mesh [1 x ndim]
+    padding_distance: list
+        Padding distances [[W,E], [N,S], [Down,Up]]
+    base_mesh: discretize.BaseMesh
+        discretize mesh used to center the core mesh
+    depth_core: float
+        Depth of core mesh below xyz
+    expansion_factor: float
+        Expension factor for padding cells [1.3]
+    mesh_type: str
+        Specify output mesh type: ["TENSOR"] or "TREE"
+
+    Returns
+    --------
+    discretize.BaseMesh
+        Mesh of "mesh_type"
+
+    Example
+    -------
+    .. plot::
+        :include-source:
+
+        import discretize
+        import matplotlib.pyplot as plt
+        import numpy as np
+
+        xyLoc = np.random.randn(8,2)
+
+        mesh = discretize.utils.meshutils.mesh_builder_xyz(
+            xyLoc, [0.1, 0.1], depth_core=0.5,
+            padding_distance=[[1,2], [1,0]],
+            mesh_type='tensor',
+        )
+
+
+        axs = plt.subplot()
+        mesh.plotImage(mesh.vol, grid=True, ax=axs)
+        axs.scatter(xyLoc[:,0], xyLoc[:,1], 15, c='w', zorder=3)
+        axs.set_aspect('equal')
+        plt.show()
+
+    """
+    if mesh_type.lower() not in ['tensor', 'tree']:
+        raise ValueError(
+            'Revise mesh_type. Only TENSOR | TREE mesh are implemented'
+        )
+
+    # Get extent of points
+    limits = []
+    center = []
+    nC = []
+    for dim in range(xyz.shape[1]):
+
+        max_min = np.r_[xyz[:, dim].max(), xyz[:, dim].min()]
+        limits += [max_min]
+        center += [np.mean(max_min)]
+        nC += [int((max_min[0] - max_min[1]) / h[dim])]
+
+    if depth_core is not None:
+        nC[-1] += int(depth_core / h[-1])
+        limits[-1][1] -= depth_core
+
+    if mesh_type.lower() == 'tensor':
+
+        # Figure out padding cells from distance
+        def expand(dx, pad):
+            length = 0
+            nc = 0
+            while length < pad:
+                nc += 1
+                length = np.sum(dx * expansion_factor**(np.asarray(range(nc))+1))
+
+            return nc
+
+        # Define h along each dimension
+        h_dim = []
+        nC_x0 = []
+        for dim in range(xyz.shape[1]):
+            h_dim += [[
+                (
+                    h[dim],
+                    expand(h[dim], padding_distance[dim][0]),
+                    -expansion_factor
+                ),
+                (h[dim], nC[dim]),
+                (
+                    h[dim],
+                    expand(h[dim], padding_distance[dim][1]),
+                    expansion_factor
+                )
+            ]]
+
+            nC_x0 += [h_dim[-1][0][1]]
+
+        # Create mesh
+        mesh = discretize.TensorMesh(h_dim)
+
+    elif mesh_type.lower() == 'tree':
+
+        # Figure out full extent required from input
+        h_dim = []
+        nC_x0 = []
+        for ii, cc in enumerate(nC):
+            extent = (
+                limits[ii][0] -
+                limits[ii][1] +
+                np.sum(padding_distance[ii])
+            )
+
+            # Number of cells at the small octree level
+            maxLevel = int(np.log2(extent / h[ii])) + 1
+            h_dim += [np.ones(2**maxLevel) * h[ii]]
+
+        # Define the mesh and origin
+        mesh = discretize.TreeMesh(h_dim)
+
+        for ii, cc in enumerate(nC):
+            core = limits[ii][0] - limits[ii][1]
+            pad2 = int(np.log2(padding_distance[ii][0] / h[ii] + 1))
+
+            nC_x0 += [int(np.ceil((mesh.h[ii].sum() - core) / h[ii] / 2))]
+
+
+    # Set origin
+    x0 = []
+    for ii, hi in enumerate(mesh.h):
+        x0 += [limits[ii][1] - np.sum(hi[:nC_x0[ii]])]
+
+    mesh.x0 = np.hstack(x0)
+
+    # Shift mesh if global mesh is used based on closest to centroid
+    axis = ["x", "y", "z"]
+    if base_mesh is not None:
+        for dim in range(base_mesh.dim):
+
+            cc_base = getattr(
+                    base_mesh,
+                    "vectorCC{orientation}".format(
+                            orientation=axis[dim]
+                        )
+            )
+
+            cc_local = getattr(
+                    mesh,
+                    "vectorCC{orientation}".format(
+                            orientation=axis[dim]
+                        )
+            )
+
+            shift = (
+                cc_base[np.max([np.searchsorted(cc_base, center[dim])-1, 0])] -
+                cc_local[np.max([np.searchsorted(cc_local, center[dim])-1, 0])]
+            )
+
+            x0[dim] += shift
+
+            mesh.x0 = np.hstack(x0)
+
+    return mesh
+
+
+def refine_tree_xyz(
+    mesh, xyz,
+    method="radial",
+    octree_levels=[1, 1, 1],
+    octree_levels_padding=None,
+    finalize=False,
+    min_level=0,
+    max_distance=np.inf
+):
+    """
+    Refine a TreeMesh based on xyz point locations
+
+    Parameters
+    ----------
+    mesh: BaseMesh
+        The TreeMesh object to be refined
+    xyz: numpy.ndarray
+        2D array of points
+    method: str
+        Method used to refine the mesh based on xyz locations
+
+        - "radial": Based on radial distance xyz and cell centers
+        - "surface": Along triangulated surface repeated vertically
+        - "box": Inside limits defined by outer xyz locations
+
+    octree_levels: list
+        Minimum number of cells around points in each k octree level
+        [N(k), N(k-1), ...]
+    octree_levels_padding: list
+        Padding cells added to the outer limits of the data each octree levels
+        used for method= "surface" and "box" [N(k), N(k-1), ...].
+    finalize: bool
+        True | [False]    Finalize the TreeMesh
+    max_distance: float
+        Maximum refinement distance from xyz locations.
+        Used for method="surface" to reduce interpolation distance.
+
+    Returns
+    --------
+    discretize.TreeMesh
+        mesh
+
+    """
+
+    if octree_levels_padding is not None:
+
+        if len(octree_levels_padding) != len(octree_levels):
+            raise ValueError(
+                "'octree_levels_padding' must be the length %i" % len(octree_levels)
+            )
+
+    else:
+        octree_levels_padding = np.zeros_like(octree_levels)
+
+    # Prime the refinement against large cells
+    mesh.insert_cells(
+        xyz,
+        [mesh.max_level - np.nonzero(octree_levels)[0][0]]*xyz.shape[0],
+        finalize=False
+    )
+
+    # Trigger different refine methods
+    if method.lower() == "radial":
+
+        # Build a cKDTree for fast nearest lookup
+        tree = cKDTree(xyz)
+
+        # Compute the outer limits of each octree level
+        rMax = np.cumsum(
+            mesh.hx.min() *
+            np.asarray(octree_levels) *
+            2**np.arange(len(octree_levels))
+        )
+
+        # Radial function
+        def inBall(cell):
+            xyz = cell.center
+            r, ind = tree.query(xyz)
+
+            for ii, nC in enumerate(octree_levels):
+
+                if r < rMax[ii]:
+
+                    return mesh.max_level-ii
+
+            return min_level
+
+        mesh.refine(inBall, finalize=finalize)
+
+    elif method.lower() == 'surface':
+
+        # Compute centroid
+        centroid = np.mean(xyz, axis=0)
+
+        if mesh.dim == 2:
+            rOut = np.abs(centroid[0]-xyz).max()
+            hz = mesh.hy.min()
+        else:
+            # Largest outer point distance
+            rOut = np.linalg.norm(np.r_[
+                np.abs(centroid[0]-xyz[:, 0]).max(),
+                np.abs(centroid[1]-xyz[:, 1]).max()
+                ]
+            )
+            hz = mesh.hz.min()
+
+        # Compute maximum depth of refinement
+        zmax = np.cumsum(
+            hz *
+            np.asarray(octree_levels) *
+            2**np.arange(len(octree_levels))
+        )
+
+        # Compute maximum horizontal padding offset
+        padWidth = np.cumsum(
+            mesh.hx.min() *
+            np.asarray(octree_levels_padding) *
+            2**np.arange(len(octree_levels_padding))
+        )
+
+        # Increment the vertical offset
+        zOffset = 0
+        xyPad = -1
+        depth = zmax[-1]
+        # Cycle through the Tree levels backward
+        for ii in range(len(octree_levels)-1, -1, -1):
+
+            dx = mesh.hx.min() * 2**ii
+
+            if mesh.dim == 3:
+                dy = mesh.hy.min() * 2**ii
+                dz = mesh.hz.min() * 2**ii
+            else:
+                dz = mesh.hy.min() * 2**ii
+
+            # Increase the horizontal extent of the surface
+            if xyPad != padWidth[ii]:
+                xyPad = padWidth[ii]
+
+                # Calculate expansion for padding XY cells
+                expansion_factor = (rOut + xyPad) / rOut
+                xLoc = (xyz - centroid)*expansion_factor + centroid
+
+                if mesh.dim == 3:
+                    # Create a new triangulated surface
+                    tri2D = Delaunay(xLoc[:, :2])
+                    F = interpolate.LinearNDInterpolator(tri2D, xLoc[:, 2])
+                else:
+                    F = interpolate.interp1d(xLoc[:, 0], xLoc[:, 1], fill_value='extrapolate')
+
+            limx = np.r_[xLoc[:, 0].max(), xLoc[:, 0].min()]
+            nCx = int(np.ceil((limx[0]-limx[1]) / dx))
+
+            if mesh.dim == 3:
+                limy = np.r_[xLoc[:, 1].max(), xLoc[:, 1].min()]
+                nCy = int(np.ceil((limy[0]-limy[1]) / dy))
+
+                # Create a grid at the octree level in xy
+                CCx, CCy = np.meshgrid(
+                    np.linspace(
+                        limx[1], limx[0], nCx
+                        ),
+                    np.linspace(
+                        limy[1], limy[0], nCy
+                        )
+                )
+
+                xy = np.c_[CCx.reshape(-1), CCy.reshape(-1)]
+
+                # Only keep points within triangulation
+                indexTri = tri2D.find_simplex(xy)
+
+            else:
+                xy = np.linspace(
+                        limx[1], limx[0], nCx
+                        )
+                indexTri = np.ones_like(xy, dtype='bool')
+
+            # Interpolate the elevation linearly
+            z = F(xy[indexTri != -1])
+
+            newLoc = np.c_[xy[indexTri != -1], z]
+
+            # Only keep points within max_distance
+            tree = cKDTree(xyz)
+            r, ind = tree.query(newLoc)
+
+            # Apply vertical padding for current octree level
+            zOffset = 0
+            while zOffset < depth:
+                indIn = r < (max_distance + padWidth[ii])
+                nnz = int(np.sum(indIn))
+                if nnz > 0:
+                    mesh.insert_cells(
+                        np.c_[
+                                newLoc[indIn, :2],
+                                newLoc[indIn, 2]-zOffset],
+                        np.ones(nnz)*mesh.max_level-ii,
+                        finalize=False
+                    )
+
+                zOffset += dz
+
+            depth -= dz * octree_levels[ii]
+
+        if finalize:
+            mesh.finalize()
+
+    elif method.lower() == 'box':
+
+        # Define the data extend [bottom SW, top NE]
+        bsw = np.min(xyz, axis=0)
+        tne = np.max(xyz, axis=0)
+
+        hx = mesh.hx.min()
+
+        if mesh.dim == 2:
+            hz = mesh.hy.min()
+        else:
+            hz = mesh.hz.min()
+
+        # Pre-calculate max depth of each level
+        zmax = np.cumsum(
+            hz * np.asarray(octree_levels) *
+            2**np.arange(len(octree_levels))
+        )
+
+        if mesh.dim == 2:
+            # Pre-calculate outer extent of each level
+            padWidth = np.cumsum(
+                    mesh.hx.min() *
+                    np.asarray(octree_levels_padding) *
+                    2**np.arange(len(octree_levels_padding))
+                )
+
+            # Make a list of outer limits
+            BSW = [
+                bsw - np.r_[padWidth[ii], zmax[ii]]
+                for ii, (octZ, octXY) in enumerate(
+                        zip(octree_levels, octree_levels_padding)
+                )
+            ]
+
+            TNE = [
+                tne + np.r_[padWidth[ii], zmax[ii]]
+                for ii, (octZ, octXY) in enumerate(
+                    zip(octree_levels, octree_levels_padding)
+                )
+            ]
+
+        else:
+            hy = mesh.hy.min()
+
+            # Pre-calculate outer X extent of each level
+            padWidth_x = np.cumsum(
+                    hx * np.asarray(octree_levels_padding) *
+                    2**np.arange(len(octree_levels_padding))
+                )
+
+            # Pre-calculate outer Y extent of each level
+            padWidth_y = np.cumsum(
+                    hy * np.asarray(octree_levels_padding) *
+                    2**np.arange(len(octree_levels_padding))
+                )
+
+            # Make a list of outer limits
+            BSW = [
+                bsw - np.r_[padWidth_x[ii], padWidth_y[ii], zmax[ii]]
+                for ii, (octZ, octXY) in enumerate(
+                        zip(octree_levels, octree_levels_padding)
+                )
+            ]
+
+            TNE = [
+                tne + np.r_[padWidth_x[ii], padWidth_y[ii], zmax[ii]]
+                for ii, (octZ, octXY) in enumerate(
+                    zip(octree_levels, octree_levels_padding)
+                )
+            ]
+
+        def inBox(cell):
+
+            xyz = cell.center
+
+            for ii, (nC, bsw, tne) in enumerate(zip(octree_levels, BSW, TNE)):
+
+                if np.all([xyz > bsw, xyz < tne]):
+                    return mesh.max_level-ii
+
+            return cell._level
+
+        mesh.refine(inBox, finalize=finalize)
+
+    else:
+        raise NotImplementedError(
+            "Only method= 'radial', 'surface'"
+            " or 'box' have been implemented"
+        )
+
+    return mesh
diff --git a/docs/api/index.rst b/docs/api/index.rst
index c2927ea53..aaf6717d4 100644
--- a/docs/api/index.rst
+++ b/docs/api/index.rst
@@ -101,6 +101,8 @@ Mesh Utilities
     utils.closestPoints
     utils.ExtractCoreMesh
     utils.random_model
+    utils.mesh_builder_xyz
+    utils.refine_tree_xyz
 
 
 Matrix Utilities
diff --git a/docs/conf.py b/docs/conf.py
index a9aba242c..b06b51de5 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -39,7 +39,7 @@
     'sphinx.ext.extlinks',
     'sphinx.ext.intersphinx',
     'sphinx.ext.mathjax',
-    # 'sphinx.ext.napoleon',
+    'sphinx.ext.napoleon',
     'sphinx.ext.todo',
     'sphinx.ext.viewcode',
     'matplotlib.sphinxext.plot_directive',
@@ -53,16 +53,16 @@
 
 numpydoc_class_members_toctree = False
 
-# napoleon_numpy_docstring = True
-# napoleon_include_init_with_doc = False
-# napoleon_include_private_with_doc = False
-# napoleon_include_special_with_doc = True
-# napoleon_use_admonition_for_examples = False
-# napoleon_use_admonition_for_notes = False
-# napoleon_use_admonition_for_references = False
-# napoleon_use_ivar = False
-# napoleon_use_param = True
-# napoleon_use_rtype = True
+napoleon_numpy_docstring = True
+napoleon_include_init_with_doc = False
+napoleon_include_private_with_doc = False
+napoleon_include_special_with_doc = True
+napoleon_use_admonition_for_examples = False
+napoleon_use_admonition_for_notes = False
+napoleon_use_admonition_for_references = False
+napoleon_use_ivar = False
+napoleon_use_param = True
+napoleon_use_rtype = True
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
@@ -155,7 +155,7 @@
     html_theme = 'sphinx_rtd_theme'
     html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
     pass
-except Exception, e:
+except Exception:
     html_theme = 'default'
 
 # Theme options are theme-specific and customize the look and feel of a theme
diff --git a/docs/examples/vtki_demo.py b/docs/examples/vtki_demo.py
index 95d1acbf4..f3402493e 100644
--- a/docs/examples/vtki_demo.py
+++ b/docs/examples/vtki_demo.py
@@ -8,7 +8,7 @@
 To run this example, you will need to `install vtki`_.
 
 .. _vtki: http://www.vtki.org
-.. _install vtki: http://www.vtki.org/en/latest/getting-started/installation.html
+.. _install vtki: http://docs.vtki.org/getting-started/installation.html
 
 - contributed by `@banesullivan <https://github.com/banesullivan>`_
 
diff --git a/tests/tree/test_tree_utils.py b/tests/tree/test_tree_utils.py
new file mode 100644
index 000000000..af28f9e44
--- /dev/null
+++ b/tests/tree/test_tree_utils.py
@@ -0,0 +1,136 @@
+from __future__ import print_function
+import numpy as np
+import unittest
+import discretize
+from discretize.utils import meshutils
+
+TOL = 1e-8
+np.random.seed(12)
+
+
+class TestRefineOcTree(unittest.TestCase):
+
+    def test_radial(self):
+        dx = 0.25
+        rad = 10
+        mesh = meshutils.mesh_builder_xyz(
+            np.c_[0.01, 0.01, 0.01], [dx, dx, dx],
+            depth_core=0,
+            padding_distance=[[0, 20], [0, 20], [0, 20]],
+            mesh_type='TREE',
+        )
+
+        radCell = int(np.ceil(rad/dx))
+        mesh = meshutils.refine_tree_xyz(
+            mesh, np.c_[0, 0, 0],
+            octree_levels=[radCell], method='radial',
+            finalize=True
+        )
+
+        # Volume of sphere
+        vol = 4.*np.pi/3.*rad**3.
+
+        residual = np.abs(
+            vol -
+            mesh.vol[
+                mesh._cell_levels_by_indexes(range(mesh.nC)) == mesh.max_level
+            ].sum()
+        )/vol * 100
+
+        self.assertTrue(residual < 3)
+
+    def test_box(self):
+        dx = 0.25
+        dl = 10
+
+        # Create a box 2*dl in width
+        X, Y, Z = np.meshgrid(np.c_[-dl, dl], np.c_[-dl, dl], np.c_[-dl, dl])
+        xyz = np.c_[np.ravel(X), np.ravel(Y), np.ravel(Z)]
+        mesh = meshutils.mesh_builder_xyz(
+            np.c_[0.01, 0.01, 0.01], [dx, dx, dx],
+            depth_core=0,
+            padding_distance=[[0, 20], [0, 20], [0, 20]],
+            mesh_type='TREE',
+        )
+
+        mesh = meshutils.refine_tree_xyz(
+            mesh, xyz, octree_levels=[1], method='box', finalize=True
+        )
+
+        # Volume of box
+        vol = (2*dl)**3
+
+        residual = np.abs(
+            vol -
+            mesh.vol[
+                mesh._cell_levels_by_indexes(range(mesh.nC)) == mesh.max_level
+            ].sum()
+        )/vol * 100
+
+        self.assertTrue(residual < 0.5)
+
+    def test_surface(self):
+        dx = 0.1
+        dl = 20
+
+        # Define triangle
+        xyz = np.r_[
+            np.c_[-dl/2, -dl/2, 0],
+            np.c_[dl/2, -dl/2, 0],
+            np.c_[dl/2, dl/2, 0]
+        ]
+        mesh = meshutils.mesh_builder_xyz(
+            np.c_[0.01, 0.01, 0.01], [dx, dx, dx],
+            depth_core=0,
+            padding_distance=[[0, 20], [0, 20], [0, 20]],
+            mesh_type='TREE',
+        )
+
+        mesh = meshutils.refine_tree_xyz(
+            mesh, xyz, octree_levels=[1], method='surface', finalize=True
+        )
+
+        # Volume of triangle
+        vol = dl*dl*dx/2
+
+        residual = np.abs(
+            vol -
+            mesh.vol[
+                mesh._cell_levels_by_indexes(range(mesh.nC)) == mesh.max_level
+            ].sum()/2
+        )/vol * 100
+
+        self.assertTrue(residual < 5)
+
+    def test_errors(self):
+        dx = 0.25
+        rad = 10
+        self.assertRaises(ValueError, meshutils.mesh_builder_xyz,
+            np.c_[0.01, 0.01, 0.01], [dx, dx, dx],
+            depth_core=0,
+            padding_distance=[[0, 20], [0, 20], [0, 20]],
+            mesh_type='cyl',
+        )
+
+        mesh = meshutils.mesh_builder_xyz(
+            np.c_[0.01, 0.01, 0.01], [dx, dx, dx],
+            depth_core=0,
+            padding_distance=[[0, 20], [0, 20], [0, 20]],
+            mesh_type='tree',
+        )
+
+        radCell = int(np.ceil(rad/dx))
+        self.assertRaises(NotImplementedError, meshutils.refine_tree_xyz,
+            mesh, np.c_[0, 0, 0], octree_levels=[radCell], method='other',
+            finalize=True
+        )
+
+        self.assertRaises(ValueError, meshutils.refine_tree_xyz,
+            mesh, np.c_[0, 0, 0], octree_levels=[radCell],
+            octree_levels_padding=[], method='surface',
+            finalize=True
+        )
+
+
+if __name__ == '__main__':
+    unittest.main()