diff --git a/.gitignore b/.gitignore
index 74f6527809..31a73b938a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 *~
 *.py~
 *.swp
+*.py#
diff --git a/src/porepy/utils/__init__.py b/src/porepy/utils/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/src/porepy/utils/accumarray.py b/src/porepy/utils/accumarray.py
new file mode 100644
index 0000000000..84bc1caa88
--- /dev/null
+++ b/src/porepy/utils/accumarray.py
@@ -0,0 +1,108 @@
+from itertools import product
+import numpy as np
+
+
+def accum(accmap, a, func=None, size=None, fill_value=0, dtype=None):
+    """
+    An accumulation function similar to Matlab's `accumarray` function.
+
+    Parameters
+    ----------
+    accmap : ndarray
+        This is the "accumulation map".  It maps input (i.e. indices into
+        `a`) to their destination in the output array.  The first `a.ndim`
+        dimensions of `accmap` must be the same as `a.shape`.  That is,
+        `accmap.shape[:a.ndim]` must equal `a.shape`.  For example, if `a`
+        has shape (15,4), then `accmap.shape[:2]` must equal (15,4).  In this
+        case `accmap[i,j]` gives the index into the output array where
+        element (i,j) of `a` is to be accumulated.  If the output is, say,
+        a 2D, then `accmap` must have shape (15,4,2).  The value in the
+        last dimension give indices into the output array. If the output is
+        1D, then the shape of `accmap` can be either (15,4) or (15,4,1)
+    a : ndarray
+        The input data to be accumulated.
+    func : callable or None
+        The accumulation function.  The function will be passed a list
+        of values from `a` to be accumulated.
+        If None, numpy.sum is assumed.
+    size : ndarray or None
+        The size of the output array.  If None, the size will be determined
+        from `accmap`.
+    fill_value : scalar
+        The default value for elements of the output array.
+    dtype : numpy data type, or None
+        The data type of the output array.  If None, the data type of
+        `a` is used.
+
+    Returns
+    -------
+    out : ndarray
+        The accumulated results.
+
+        The shape of `out` is `size` if `size` is given.  Otherwise the
+        shape is determined by the (lexicographically) largest indices of
+        the output found in `accmap`.
+
+
+    Examples
+    --------
+    >>> from numpy import array, prod
+    >>> a = array([[1,2,3],[4,-1,6],[-1,8,9]])
+    >>> a
+    array([[ 1,  2,  3],
+           [ 4, -1,  6],
+           [-1,  8,  9]])
+    >>> # Sum the diagonals.
+    >>> accmap = array([[0,1,2],[2,0,1],[1,2,0]])
+    >>> s = accum(accmap, a)
+    array([9, 7, 15])
+    >>> # A 2D output, from sub-arrays with shapes and positions like this:
+    >>> # [ (2,2) (2,1)]
+    >>> # [ (1,2) (1,1)]
+    >>> accmap = array([
+            [[0,0],[0,0],[0,1]],
+            [[0,0],[0,0],[0,1]],
+            [[1,0],[1,0],[1,1]],
+        ])
+    >>> # Accumulate using a product.
+    >>> accum(accmap, a, func=prod, dtype=float)
+    array([[ -8.,  18.],
+           [ -8.,   9.]])
+    >>> # Same accmap, but create an array of lists of values.
+    >>> accum(accmap, a, func=lambda x: x, dtype='O')
+    array([[[1, 2, 4, -1], [3, 6]],
+           [[-1, 8], [9]]], dtype=object)
+    """
+
+    # Check for bad arguments and handle the defaults.
+    if accmap.shape[:a.ndim] != a.shape:
+        raise ValueError("The initial dimensions of accmap must be the same as a.shape")
+    if func is None:
+        func = np.sum
+    if dtype is None:
+        dtype = a.dtype
+    if accmap.shape == a.shape:
+        accmap = np.expand_dims(accmap, -1)
+    adims = tuple(range(a.ndim))
+    if size is None:
+        size = 1 + np.squeeze(np.apply_over_axes(np.max, accmap, axes=adims))
+    size = np.atleast_1d(size)
+
+    # Create an array of python lists of values.
+    vals = np.empty(size, dtype='O')
+    for s in product(*[range(k) for k in size]):
+        vals[s] = []
+    for s in product(*[range(k) for k in a.shape]):
+        indx = tuple(accmap[s])
+        val = a[s]
+        vals[indx].append(val)
+
+    # Create the output array.
+    out = np.empty(size, dtype=dtype)
+    for s in product(*[range(k) for k in size]):
+        if vals[s] == []:
+            out[s] = fill_value
+        else:
+            out[s] = func(vals[s])
+
+    return out
diff --git a/src/porepy/utils/errors/error.py b/src/porepy/utils/errors/error.py
new file mode 100644
index 0000000000..df8b5d8540
--- /dev/null
+++ b/src/porepy/utils/errors/error.py
@@ -0,0 +1,106 @@
+import numpy as np
+
+#------------------------------------------------------------------------------#
+
+def interpolate(g, fun):
+    """
+    Interpolate a scalar or vector function on the cell centers of the grid.
+
+    Parameters
+    ----------
+    g : grid
+        Grid, or a subclass, with geometry fields computed.
+    fun : function
+        Scalar or vector function.
+
+    Return
+    ------
+    out: np.ndarray (dim of fun, g.num_cells)
+        Function interpolated in the cell centers.
+
+    Examples
+    --------
+    def fun_p(pt): return np.sin(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1])
+    def fun_u(pt): return [\
+                      -2*np.pi*np.cos(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1]),
+                      -2*np.pi*np.sin(2*np.pi*pt[0])*np.cos(2*np.pi*pt[1])]
+    p_ex = interpolate(g, fun_p)
+    u_ex = interpolate(g, fun_u)
+
+    """
+
+    return np.array([fun(pt) for pt in g.cell_centers.T]).T
+
+#------------------------------------------------------------------------------#
+
+def norm_L2(g, val):
+    """
+    Compute the L2 norm of a scalar or vector field.
+
+    Parameters
+    ----------
+    g : grid
+        Grid, or a subclass, with geometry fields computed.
+    val : np.ndarray (dim of val, g.num_cells)
+        Scalar or vector field.
+
+    Return
+    ------
+    out: double
+        The L2 norm of the input field.
+
+    Examples
+    --------
+    def fun_p(pt): return np.sin(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1])
+    p_ex = interpolate(g, fun_p)
+    norm_ex = norm_L2(g, p_ex)
+
+    """
+
+    val = np.asarray(val)
+    norm_sq = lambda v: np.sum(np.multiply(np.square(v), g.cell_volumes))
+    if val.ndim == 1:
+        return np.sqrt(norm_sq(val))
+    else:
+        return np.sqrt(np.sum([norm_sq(v) for v in val]))
+
+#------------------------------------------------------------------------------#
+
+def error_L2(g, val, val_ex, relative=True):
+    """
+    Compute the L2 error of a scalar or vector field with respect to a reference
+    field. It is possible to compute the relative error (default) or the
+    absolute error.
+
+    Parameters
+    ----------
+    g : grid
+        Grid, or a subclass, with geometry fields computed.
+    val : np.ndarray (dim of val, g.num_cells)
+        Scalar or vector field.
+    val_ex : np.ndarray (dim of val, g.num_cells)
+        Reference scalar or vector field, i.e. the exact solution
+    relative: bool (True default)
+        Compute the relative error (if True) or the absolute error (if False).
+
+    Return
+    ------
+    out: double
+        The L2 error of the input fields.
+
+    Examples
+    --------
+    p = ...
+    def fun_p(pt): return np.sin(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1])
+    p_ex = interpolate(g, fun_p)
+    err_p = err_L2(g, p, p_ex)
+
+    """
+
+    val, val_ex = np.asarray(val), np.asarray(val_ex)
+    err = norm_L2(g, np.subtract(val, val_ex))
+    den = norm_L2(g, val_ex) if relative else 1
+    assert den != 0
+    return err/den
+
+#------------------------------------------------------------------------------#
diff --git a/src/porepy/utils/graph.py b/src/porepy/utils/graph.py
new file mode 100644
index 0000000000..89fb34b4b3
--- /dev/null
+++ b/src/porepy/utils/graph.py
@@ -0,0 +1,59 @@
+import numpy as np
+from scipy import sparse as sps
+
+
+class Graph:
+    """
+    A graph class.
+
+    The graph class stores the nodes and edges of the graph in a sparse
+    array (equivalently to face_nodes in the Grid class).
+
+    Attributes:
+        node_connections (sps.csc-matrix): Should be given at construction.
+            node_node connections. Matrix size: num_nodes x num_nodes.
+            node_connections[i,j] should be true if there is an edge
+            connecting node i and j.
+        regions (int) the number of regions. A region is a set of nodes
+            that can be reached by traversing the graph. Two nodes are
+            int different regions if they can not be reached by traversing
+            the graph.
+        color (int) the color of each region. Initialized as (NaN). By
+            calling color_nodes() all nodes in a region are given the
+            same colors and nodes in different regions are given different
+            colors.
+    """
+
+    def __init__(self, node_connections):
+        self.node_connections = node_connections
+        self.regions = 0
+        self.color = np.array([np.NaN] * node_connections.shape[1])
+
+    def color_nodes(self):
+        """
+        Color the nodes in each region
+        """
+        color = 0
+        while self.regions <= self.node_connections.shape[1]:
+            start = np.ravel(np.argwhere(np.isnan(self.color)))
+            if start.size != 0:
+                self.bfs(start[0], color)
+                color += 1
+                self.regions += 1
+            else:
+                return
+        raise RuntimeWarning('number of regions can not be greater than '
+                             'number of nodes')
+
+    def bfs(self, start, color):
+        """
+        Breadth first search
+        """
+        visited, queue = [], [start]
+        while queue:
+            node = queue.pop(0)
+            if node not in visited:
+                visited.append(node)
+                (_, neighbours, _) = sps.find(self.node_connections[node, :])
+                queue.extend(neighbours)
+        self.color[visited] = color
diff --git a/src/porepy/utils/half_space.py b/src/porepy/utils/half_space.py
new file mode 100644
index 0000000000..1ead73475a
--- /dev/null
+++ b/src/porepy/utils/half_space.py
@@ -0,0 +1,104 @@
+import numpy as np
+import scipy.optimize as opt
+
+
+def half_space_int(n,x0, pts):
+    """
+    Find the points that lie in the intersection of half spaces (3D)
+
+    Parameters
+    ----------
+    n : ndarray
+        This is the normal vectors of the half planes. The normal
+        vectors is assumed to point out of the half spaces.
+    x0 : ndarray
+        Point on the boundary of the half-spaces. Half space i is given
+        by all points satisfying (x - x0[:,i])*n[:,i]<=0
+    pts : ndarray
+        The points to be tested if they are in the intersection of all
+        half-spaces or not.
+
+    Returns
+    -------
+    out : ndarray
+        A logical array with length equal number of pts.
+
+        out[i] is True if pts[:,i] is in all half-spaces
+
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> n = np.array([[0,1],[1,0],[0,0]])
+    >>> x0 = np.array([[0,-1],[0,0],[0,0]])
+    >>> pts = np.array([[-1,-1,4],[2,-2,-2],[0,0,0]])
+    >>> half_space_int(n,x0,pts)
+    array([False,  True, False], dtype=bool)
+    """
+    assert n.shape[0] == 3, ' only 3D supported'
+    assert x0.shape[0] == 3, ' only 3D supported'
+    assert pts.shape[0] == 3, ' only 3D supported'
+    assert n.shape[1] == x0.shape[1], 'ther must be the same number of normal vectors as points'
+
+    n_pts   = pts.shape[1]
+    in_hull = np.zeros(n_pts)
+    x0      = np.repeat(x0[:,:,np.newaxis], n_pts,axis=2)
+    n       = np.repeat( n[:,:,np.newaxis], n_pts,axis=2)
+    for i in range(x0.shape[1]):
+        in_hull += np.sum((pts - x0[:,i,:])*n[:,i,:],axis=0)<=0
+
+    return in_hull==x0.shape[1]
+
+#------------------------------------------------------------------------------#
+
+def half_space_pt(n, x0, pts, recompute=True):
+    """
+    Find an interior point for the halfspaces.
+
+    Parameters
+    ----------
+    n : ndarray
+        This is the normal vectors of the half planes. The normal
+        vectors is assumed to coherently for all the half spaces
+        (inward or outward).
+    x0 : ndarray
+        Point on the boundary of the half-spaces. Half space i is given
+        by all points satisfying (x - x0[:,i])*n[:,i]<=0
+    pts : ndarray
+        Points which defines a bounds for the algorithm.
+    recompute: bool
+        If the algorithm fails try again with flipped normals.
+
+    Returns
+    -------
+    out: array
+        Interior point of the halfspaces.
+
+    We use linear programming to find one interior point for the half spaces.
+    Assume, n halfspaces defined by: aj*x1+bj*x2+cj*x3+dj<=0, j=1..n.
+    Perform the following linear program:
+    max(x5) aj*x1+bj*x2+cj*x3+dj*x4+x5<=0, j=1..n
+
+    Then, if [x1,x2,x3,x4,x5] is an optimal solution with x4>0 and x5>0 we get:
+    aj*(x1/x4)+bj*(x2/x4)+cj*(x3/x4)+dj<=(-x5/x4) j=1..n and (-x5/x4)<0,
+    and conclude that the point [x1/x4,x2/x4,x3/x4] is in the interior of all
+    the halfspaces. Since x5 is optimal, this point is "way in" the interior
+    (good for precision errors).
+    http://www.qhull.org/html/qhalf.htm#notes
+
+    """
+    dim = (1,n.shape[1])
+    c = np.array([0,0,0,0,-1])
+    A_ub = np.concatenate( (n, [np.sum(-n*x0, axis=0)], np.ones(dim)) ).T
+    bounds = ( (np.amin(pts[0,:]), np.amax(pts[0,:]) ),
+               (np.amin(pts[1,:]), np.amax(pts[1,:]) ),
+               (np.amin(pts[2,:]), np.amax(pts[2,:]) ),
+               (None, None), (None, None) )
+    res = opt.linprog(c, A_ub, np.zeros(dim).T, bounds=bounds)
+    if recompute and (res.status != 0 or res.x[3] <=0 or res.x[4] <= 0):
+        return half_space_pt(-n, x0, pts, False)
+    else:
+        assert res.status == 0 and res.x[3] > 0 and res.x[4] > 0
+        return np.array(res.x[0:3]) / res.x[3]
+
+#------------------------------------------------------------------------------#
diff --git a/src/porepy/utils/matrix_compression.py b/src/porepy/utils/matrix_compression.py
new file mode 100644
index 0000000000..ff28641a2c
--- /dev/null
+++ b/src/porepy/utils/matrix_compression.py
@@ -0,0 +1,50 @@
+""" Functions for compressing matrices to compact format, and recover them.
+
+Acknowledgements:
+    The functions are a python translation of the corresponding matlab
+    functions found in the Matlab Reservoir Simulation Toolbox (MRST) developed
+    by SINTEF ICT, see www.sintef.no/projectweb/mrst/ . 
+    
+"""
+
+import numpy as np
+
+
+def rldecode(A,n):
+    """ Decode compressed information. 
+        
+        The code is heavily inspired by MRST's function with the same name, 
+        however, requirements on the shape of functions are probably somewhat
+        different.
+        
+        >>> rldecode(np.array([1, 2, 3]), np.array([2, 3, 1]))
+        [1, 1, 2, 2, 2, 3]
+        
+        >>> rldecode(np.array([1, 2]), np.array([1, 3]))
+        [1, 2, 2, 2]
+        
+        Args:
+            A (double, m x k), compressed matrix to be recovered. The 
+            compression should be along dimension 1
+            n (int): Number of occurences for each element
+    """
+    r = n > 0         
+    i = np.cumsum(np.hstack((np.zeros(1), n[r])), dtype='>i4')
+    j = np.zeros(i[-1])
+    j[i[1:-1:]] = 1
+    B = A[np.cumsum(j, dtype='>i4')]
+    return B
+
+def rlencode(A):
+    """ Compress matrix by looking for identical columns. """
+    comp = A[::, 0:-1] != A[::, 1::]
+    i = np.any(comp, axis=0)
+    i = np.hstack((np.argwhere(i).ravel(), (A.shape[1]-1)))
+    
+    num = np.diff(np.hstack((np.array([-1]), i)))
+    
+    return A[::, i], num
+    
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()
diff --git a/src/porepy/utils/mcolon.py b/src/porepy/utils/mcolon.py
new file mode 100644
index 0000000000..9ef700d632
--- /dev/null
+++ b/src/porepy/utils/mcolon.py
@@ -0,0 +1,65 @@
+""" Efficient numpy.arange for arrays of start and end indices.
+
+Acknowledgements:
+    The functions are a python translation of the corresponding matlab
+    functions found in the Matlab Reservoir Simulation Toolbox (MRST) developed
+    by SINTEF ICT, see www.sintef.no/projectweb/mrst/ .
+
+"""
+
+import numpy as np
+
+
+def mcolon(lo, hi):
+    """ Expansion of np.arange(a, b) for arrays a and b.
+
+    The code is equivalent to the following (less efficient) loop:
+    arr = np.empty(0)
+    for l, h in zip(lo, hi):
+        arr = np.hstack((arr, np.arange(l, h+1, 1)))
+
+    Parameters:
+        lo (np.ndarray, int): Lower bounds of the arrays to be created.
+        hi (np.ndarray, int): Upper bounds of the arrays to be created. The
+            elements in hi will be included in the resulting array.
+
+        lo and hi should either have 1 or n elements. If their size are both
+        larger than one, they should have the same length.
+
+    Examples:
+        >>> mcolon(np.array([0, 0, 0]), np.array([1, 3, 2]))
+        array([0, 1, 0, 1, 2, 3, 0, 1, 2], dtype=int64)
+
+        >>> mcolon(np.array([0, 1]), np.array([2]))
+        array([0, 1, 2, 1, 2], dtype=int64)
+
+        >>> mcolon(np.array([0, 1, 1, 1]), np.array([0, 3, 3, 3]))
+        array([0, 1, 2, 3, 1, 2, 3, 1, 2, 3], dtype=int64)
+
+    """
+    if lo.size == 1:
+        lo = lo * np.ones(hi.size, dtype='int64')
+    if hi.size == 1:
+        hi = hi * np.ones(lo.size, dtype='int64')
+    if lo.size != hi.size:
+        raise ValueError('Low and high should have same number of elements, '
+                         'or a single item ')
+
+    i = hi >= lo
+    if not any(i):
+        return None
+
+    lo = lo[i]
+    hi = hi[i]
+    d = hi - lo + 1
+    n = np.sum(d)
+
+    x = np.ones(n, dtype='int64')
+    x[0] = lo[0]
+    x[np.cumsum(d[0:-1]).astype('int64')] = lo[1:] - hi[0:-1]
+    return np.cumsum(x)
+
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()
diff --git a/src/porepy/utils/permutations.py b/src/porepy/utils/permutations.py
new file mode 100644
index 0000000000..48aec5b164
--- /dev/null
+++ b/src/porepy/utils/permutations.py
@@ -0,0 +1,46 @@
+
+def multinary_permutations(base, length):
+    """
+    Define a generator over all numbers of a certain length for a number system
+    with a specified base.
+
+    For details on the decomposition into an arbitrary base see
+        http://math.stackexchange.com/questions/111150/changing-a-number-between-arbitrary-bases
+
+    Note that the generator will loop over base**length combinations.
+
+    Examples:
+
+        Construct the numbers [00] to [11] in binary numbers
+        >>> multinary_permutations(2, 2)
+        [array([ 0.,  0.]), array([ 1.,  0.]), array([ 0.,  1.]), array([ 1.,  1.])]
+
+        Construct the numbers from 0 to 99 (permuted) in the decimal number
+        system.
+        >>> it = multinary_permutation(10, 2)
+
+    Parameters:
+        base (int): Base of the number system
+        length (int): Number of digits in the numbers
+
+    Yields:
+        array, size length: Array describing the next number combination.
+
+    """
+
+    # There are in total base ** length numbers to be covered, these need to be
+    # rewritten into the base number system
+    for iter1 in range(base **  length):
+
+        # Array to store the multi-d index of the current index
+        bit_val = [0] * length
+        # Number to be decomposed
+        v = iter1
+
+        # Loop over all digits, find the expression of v in that system
+        for iter2 in range(length):
+            bit_val[iter2] = v % base
+            v = v // base
+        # Yield the next value
+        yield bit_val
+
diff --git a/src/porepy/utils/setmembership.py b/src/porepy/utils/setmembership.py
new file mode 100644
index 0000000000..f3debd3f76
--- /dev/null
+++ b/src/porepy/utils/setmembership.py
@@ -0,0 +1,278 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 25 20:16:42 2016
+
+@author: keile
+"""
+
+import numpy as np
+
+import utils.unique as np_unique
+
+def unique_rows(data):
+    """
+    Function similar to Matlab's unique(...,'rows')
+
+    See also function unique_columns in this module; this is likely slower, but
+    is understandable, documented, and has a tolerance option.
+
+    Copied from
+    http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array/
+    (summary pretty far down on the page)
+    Note: I have no idea what happens here
+
+    """
+    b = np.ascontiguousarray(data).view(np.dtype((np.void,
+                                                  data.dtype.itemsize * data.shape[1])))
+    _, ia = np.unique(b, return_index=True)
+    _, ic = np.unique(b, return_inverse=True)
+#    return np.unique(b).view(data.dtype).reshape(-1, data.shape[1]), ia, ic
+    return data[ia], ia, ic
+
+
+def _asvoid(arr):
+    """
+
+    Taken from
+    http://stackoverflow.com/questions/22699756/python-version-of-ismember-with-rows-and-index
+
+    View the array as dtype np.void (bytes)
+    This views the last axis of ND-arrays as bytes so you can perform
+    comparisons on the entire row.
+    http://stackoverflow.com/a/16840350/190597 (Jaime, 2013-05)
+    Warning: When using asvoid for comparison, note that float zeros may
+    compare UNEQUALLY
+    >>> asvoid([-0.]) == asvoid([0.])
+    array([False], dtype=bool)
+    """
+    arr = np.ascontiguousarray(arr)
+    return arr.view(np.dtype((np.void, arr.dtype.itemsize * arr.shape[-1])))
+
+
+def _find_occ(a, b):
+    """
+    Find index of occurences of a in b.
+
+    The function has only been tested on np.arrays, but it should be fairly
+    general (only require iterables?)
+
+    Code snippet found at
+    http://stackoverflow.com/questions/15864082/python-equivalent-of-matlabs-ismember-function?rq=1
+
+    """
+    # Base search on a dictionary
+
+    bind = {}
+    # Invert dictionary to create a map from an item in b to the *first*
+    # occurence of the item.
+    # NOTE: If we ever need to give the option of returning last index, it
+    # should require no more than bypassing the if statement.
+    for i, elt in enumerate(b):
+        if elt not in bind:
+            bind[elt] = i
+    # Use inverse mapping to obtain
+    return [bind.get(itm, None) for itm in a]
+
+
+def ismember_rows(a, b, sort=True, simple_version=False):
+    """
+    Find *columns* of a that are also members of *columns* of b.
+
+    The function mimics Matlab's function ismember(..., 'rows').
+
+    TODO: Rename function, this is confusing!
+
+    Parameters:
+        a (np.array): Each column in a will search for an equal in b.
+        b (np.array): Array in which we will look for a twin
+        sort (boolean, optional): If true, the arrays will be sorted before
+            seraching, increasing the chances for a match. Defaults to True.
+        simple_verion (boolean, optional): Use an alternative implementation
+            based on a global for loop. The code is slow for large arrays, but
+            easy to understand. Defaults to False.
+
+    Returns:
+        np.array (boolean): For each column in a, true if there is a
+            corresponding column in b.
+        np.array (int): Indexes so that b[:, ind] is also found in a.
+
+    Examples:
+        >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]])
+        >>> b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]])
+        >>> ismember_rows(a, b)
+        (array([ True,  True,  True,  True, False], dtype=bool), [1, 0, 2, 1])
+
+        >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]])
+        >>> b = np.array([[3, 1, 2, 5, 1], [3, 3, 3, 1, 2]])
+        >>> ismember_rows(a, b, sort=False)
+        (array([ True,  True, False,  True, False], dtype=bool), [1, 0, 1])
+
+    """
+
+    # Sort if required, but not if the input is 1d
+    if sort and a.ndim > 1:
+        sa = np.sort(a, axis=0)
+        sb = np.sort(b, axis=0)
+    else:
+        sa = a
+        sb = b
+
+    num_a = a.shape[-1]
+
+    if simple_version:
+        # Use straightforward search, based on a for loop. This is slow for
+        # large arrays, but as the alternative implementation is opaque, and
+        # there has been some doubts on its reliability, this version is kept
+        # as a safeguard.
+        ismem_a = np.zeros(num_a, dtype=np.bool)
+        ind_of_a_in_b = np.empty(0)
+        for i in range(num_a):
+            if sa.ndim == 1:
+                diff = np.abs(sb - sa[i])
+            else:
+                diff = np.sum(np.abs(sb - sa[:, i].reshape((-1, 1))), axis=0)
+            if np.any(diff == 0):
+                ismem_a[i] = True
+                hit = np.where(diff == 0)[0]
+                if hit.size > 1:
+                    hit = hit[0]
+                ind_of_a_in_b = np.append(ind_of_a_in_b, hit)
+
+        return ismem_a, ind_of_a_in_b.astype('int')
+
+    else:
+        if a.ndim == 1:
+            # Special treatment of 1d, vstack of voids (below) ran into trouble
+            # here.
+            unq, k, count = np.unique(np.hstack((a, b)), return_inverse=True,
+                                      return_counts=True)
+            _, k_a, count_a = np.unique(a, return_inverse=True,
+                                        return_counts=True)
+        else:
+            # Represent the arrays as voids to facilitate quick comparison
+            # Take void type of int64s, or else spurious error messages may
+            # arise. We do this after the transpose (which copies sa and sb) to
+            # avoid disturbing the original fields.
+            voida = _asvoid(sa.transpose().astype('int64'))
+            voidb = _asvoid(sb.transpose().astype('int64'))
+
+            # Use unique to count the number of occurences in a
+            unq, j, k, count = np.unique(np.vstack((voida, voidb)),
+                                         return_index=True,
+                                         return_inverse=True,
+                                         return_counts=True)
+            # Also count the number of occurences in voida
+            _, j_a, k_a, count_a = np.unique(voida, return_index=True,
+                                             return_inverse=True,
+                                             return_counts=True)
+
+        # Index of a and b elements in the combined array
+        ind_a = np.arange(num_a)
+        ind_b = num_a + np.arange(b.shape[-1])
+
+        # Count number of occurences in combine array, and in a only
+        num_occ_a_and_b = count[k[ind_a]]
+        num_occ_a = count_a[k_a[ind_a]]
+
+        # Subtraction gives number of a in b
+        num_occ_a_in_b = num_occ_a_and_b - num_occ_a
+        ismem_a = (num_occ_a_in_b > 0)
+
+        # To get the indices of common elements in a and b, compare the
+        # elements in k (pointers to elements in the unique combined arary)
+        occ_a = k[ind_a[ismem_a]]
+        occ_b = k[ind_b]
+
+        ind_of_a_in_b = _find_occ(occ_a, occ_b)
+        # Remove None types when no hit was found
+        ind_of_a_in_b = [i for i in ind_of_a_in_b if i is not None]
+
+        return ismem_a, np.array(ind_of_a_in_b, dtype='int')
+
+#---------------------------------------------------------
+
+
+def unique_columns_tol(mat, tol=1e-8, exponent=2):
+    """
+    Remove duplicates from a point set, for a given distance traveling.
+
+    Resembles Matlab's uniquetol function, as applied to columns. To rather
+    work at rows, use a transpose. 
+
+
+    Parameters:
+        mat (np.ndarray, nd x n_pts): Columns to be uniquified
+        tol (double, optional): Tolerance for when columns are considered equal.
+            Should be seen in connection with distance between the points in
+            the points (due to rounding errors). Defaults to 1e-8.
+        exponent (double, optional): Exponnet in norm used in distance
+            calculation. Defaults to 2.
+
+    Returns:
+        np.ndarray: Unique columns.
+        new_2_old: Index of which points that are preserved
+        old_2_new: Index of the representation of old points in the reduced
+            list.
+
+    Example (won't work as doctest):
+        >>> p_un, n2o, o2n = unique_columns(np.array([[1, 0, 1], [1, 0, 1]]))
+        >>> p_un
+        array([[1, 0], [1, 0]])
+        >>> n2o
+        array([0, 1])
+        >>> o2n
+        array([0, 1, 0])
+
+    """
+
+    # Special treatment of the case with an empty array
+    if mat.shape[1] == 0:
+         return mat
+
+    # If the matrix is integers, and the tolerance less than 1/2, we can use
+    # the new unique function that ships with numpy 1.13.
+    if issubclass(mat.dtype.type, np.integer) and tol < 0.5:
+        un_ar, new_2_old, old_2_new \
+            = np_unique.unique_np1130(mat, return_index=True,
+                                      return_inverse=True, axis=1)
+        return un_ar, new_2_old, old_2_new
+
+    def dist(p, pset):
+        " Helper function to compute distance "
+        if p.ndim == 1:
+            pt = p.reshape((-1, 1))
+        else:
+            pt = p
+
+        return np.power(np.sum(np.power(np.abs(pt - pset), exponent),
+                               axis=0), 1 / exponent)
+
+    (nd, l) = mat.shape
+
+    # By default, no columns are kept
+    keep = np.zeros(l, dtype=np.bool)
+
+    # We will however keep the first point
+    keep[0] = True
+    keep_counter = 1
+
+    # Map from old points to the unique subspace. Defaults to map to itself.
+    old_2_new = np.arange(l)
+
+    # Loop over all points, check if it is already represented in the kept list
+    for i in range(1, l):
+        proximate = np.argwhere(
+            dist(mat[:, i], mat[:, keep]) < tol * np.sqrt(nd))
+
+        if proximate.size > 0:
+            # We will not keep this point
+            old_2_new[i] = proximate[0]
+        else:
+            # We have found a new point
+            keep[i] = True
+            old_2_new[i] = keep_counter
+            keep_counter += 1
+    # Finally find which elements we kept
+    new_2_old = np.argwhere(keep).ravel()
+
+    return mat[:, keep], new_2_old, old_2_new
diff --git a/src/porepy/utils/tests/__init__.py b/src/porepy/utils/tests/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/src/porepy/utils/tests/test_graph.py b/src/porepy/utils/tests/test_graph.py
new file mode 100644
index 0000000000..48aab587f5
--- /dev/null
+++ b/src/porepy/utils/tests/test_graph.py
@@ -0,0 +1,57 @@
+import unittest
+import numpy as np
+from scipy import sparse as sps
+
+from utils import graph
+
+#------------------------------------------------------------------------------#
+
+class BasicsTest( unittest.TestCase ):
+
+#------------------------------------------------------------------------------#
+
+    def test_fully_connected_graph(self):
+        node_connections = np.array([[0,1,1,0],
+                                     [1,0,1,0],
+                                     [1,1,0,1],
+                                     [0,0,1,0]])
+        node_connections = sps.csc_matrix(node_connections)
+        G = graph.Graph(node_connections)
+        G.color_nodes()
+        assert np.all(G.color==0)
+        assert G.regions == 1
+
+#------------------------------------------------------------------------------#
+
+#------------------------------------------------------------------------------#
+
+    def test_two_region_graph(self):
+        node_connections = np.array([[0,1,1,0],
+                                     [1,0,1,0],
+                                     [1,1,0,0],
+                                     [0,0,0,0]])
+        node_connections = sps.csc_matrix(node_connections)
+        G = graph.Graph(node_connections)
+        G.color_nodes()
+        assert np.all(G.color[:2]==0)
+        assert G.color[3]==1
+        assert G.regions == 2
+
+#------------------------------------------------------------------------------#
+
+#------------------------------------------------------------------------------#
+
+    def test_nodes_connected_to_self(self):
+        node_connections = np.array([[1,1,1,0],
+                                     [1,1,1,0],
+                                     [1,1,1,0],
+                                     [0,0,0,1]])
+        node_connections = sps.csc_matrix(node_connections)
+        G = graph.Graph(node_connections)
+        G.color_nodes()
+        print(G.color)
+        assert np.all(G.color[:2]==0)
+        assert G.color[3]==1
+        assert G.regions == 2
+
+#------------------------------------------------------------------------------#
diff --git a/src/porepy/utils/tests/test_half_space.py b/src/porepy/utils/tests/test_half_space.py
new file mode 100644
index 0000000000..0e96391a9a
--- /dev/null
+++ b/src/porepy/utils/tests/test_half_space.py
@@ -0,0 +1,67 @@
+import unittest
+import numpy as np
+
+from utils import half_space
+
+#------------------------------------------------------------------------------#
+
+class BasicsTest( unittest.TestCase ):
+
+#------------------------------------------------------------------------------#
+
+    def test_one_half_space(self):
+
+        n  = np.array([[-1],[0],[0]])
+        x0 = np.array([[0],[0],[0]])
+        pts = np.array([[1,-1],[0,0],[0,0]])
+        out = half_space.half_space_int(n,x0,pts)
+        assert np.all(out == np.array([True,False]))
+
+#------------------------------------------------------------------------------#
+
+    def test_two_half_spaces(self):
+        n  = np.array([[-1,0],[0,-1],[0,0]])
+        x0 = np.array([[0,0],[0,1],[0,0]])
+        pts = np.array([[1,-1,1,0],[2,0,2,0],[0,0,0,0]])
+        out = half_space.half_space_int(n,x0,pts)
+        assert np.all(out == np.array([True,False,True,False]))
+
+#------------------------------------------------------------------------------#
+
+    def test_half_space_pt_convex_2d(self):
+        n = np.array([[0,0,0,0,1,-1],[1,1,-1,-1,0,0],[0,0,0,0,0,0]])
+        x0 = np.array([[0,0,0,0,0,2/3],[0,0,1/3,1/3,0,0],[0,0,0,0,0,0]])
+        pts = np.array([[0,2/3],[0,1/3],[0,0]])
+        pt = half_space.half_space_pt(n, x0, pts)
+        assert np.allclose(pt,[1/6,1/6,0])
+
+#------------------------------------------------------------------------------#
+
+    def test_half_space_pt_star_shaped_2d(self):
+        n = np.array([[0,0,1,0,-1,0,1],[1,1,0,1,0,-1,0],[0,0,0,0,0,0,0]])
+        x0 = np.array([[0,0,2/3,0,1,0,0],[1/3,1/3,0,0,0,2/3,0],[0,0,0,0,0,0,0]])
+        pts = np.array([[0,1],[0,2/3],[0,0]])
+        pt = half_space.half_space_pt(n, x0, pts)
+        assert np.allclose(pt,[5/6,1/2,0])
+
+#------------------------------------------------------------------------------#
+
+    def test_half_space_pt_convex_3d(self):
+        n = np.array([[1,-1,0,0,0,0],[0,0,1,-1,0,0],[0,0,0,0,1,-1]])
+        x0 = np.array([[0,1,0,0,0,0],[0,0,0,1,0,0],[0,0,0,0,0,1]])
+        pts = np.array([[1,0,0],[0,1,0],[0,0,1]])
+        pt = half_space.half_space_pt(n, x0, pts)
+        assert np.allclose(pt,[1/2,1/2,1/2])
+
+#------------------------------------------------------------------------------#
+
+    def test_half_space_pt_star_shaped_3d(self):
+        n = np.array([[0,0,1,0,-1,0,1,0,0],[1,1,0,1,0,-1,0,0,0],
+                      [0,0,0,0,0,0,0,1,-1]])
+        x0 = np.array([[0,0,2/3,0,1,0,0,0,0],[1/3,1/3,0,0,0,2/3,0,0,0],
+                       [0,0,0,0,0,0,0,0,1]])
+        pts = np.array([[0,1],[0,2/3],[0,1]])
+        pt = half_space.half_space_pt(n, x0, pts)
+        assert np.allclose(pt,[5/6,1/2,1/6])
+
+#------------------------------------------------------------------------------#
diff --git a/src/porepy/utils/tests/test_mcolon.py b/src/porepy/utils/tests/test_mcolon.py
new file mode 100644
index 0000000000..dade528d2d
--- /dev/null
+++ b/src/porepy/utils/tests/test_mcolon.py
@@ -0,0 +1,24 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Apr 14 13:51:05 2016
+
+@author: eke001
+"""
+
+import numpy as np
+
+from utils import mcolon
+
+
+def test_mcolon_simple():
+    a = np.array([1, 2])
+    b = np.array([2, 3])
+    c = mcolon.mcolon(a, b)
+    assert np.all((c - np.array([1, 2, 2, 3])) == 0)
+
+
+def test_mcolon_one_missing():
+    a = np.array([1, 2])
+    b = np.array([2, 1])
+    c = mcolon.mcolon(a, b)
+    assert np.all((c - np.array([1, 2])) == 0)
diff --git a/src/porepy/utils/tests/test_permutations.py b/src/porepy/utils/tests/test_permutations.py
new file mode 100644
index 0000000000..67359d5f46
--- /dev/null
+++ b/src/porepy/utils/tests/test_permutations.py
@@ -0,0 +1,46 @@
+import unittest
+import numpy as np
+
+from utils import permutations
+
+class TestPermutations(unittest.TestCase):
+
+    def compare_lists(self, base, length, lst):
+    # Compare a pre-defined list with the result of multinary_permutations
+    # Define a generator, and check that all values produced are contained within lst
+    # Also count the number of iterations
+        iter_cnt = 0
+        for a in permutations.multinary_permutations(base, length):
+            found = False
+            for b in lst:
+                if np.array_equal(np.array(a), np.array(b)):
+                    found = True
+                    break
+            assert found
+            iter_cnt += 1
+        assert iter_cnt == len(lst)
+
+
+    def test_length_2(self):
+        # Explicitly construct a 2D array of all combination of base 3
+        base = 3
+        length = 2
+        lst = []
+        for i in range(base):
+            for j in range(base):
+                lst.append([i, j])
+        self.compare_lists(base, length, lst)
+
+    def test_base_4(self):
+        # Check with a manually constructed list of length 3
+        base = 4
+        length = 3
+        lst = []
+        for i in range(base):
+            for j in range(base):
+                for k in range(base):
+                    lst.append([i, j, k])
+        self.compare_lists(base, length, lst)
+
+    if __name__ == '__main__':
+        unittest.main()
diff --git a/src/porepy/utils/tests/test_setmembership.py b/src/porepy/utils/tests/test_setmembership.py
new file mode 100644
index 0000000000..d63d2045d7
--- /dev/null
+++ b/src/porepy/utils/tests/test_setmembership.py
@@ -0,0 +1,169 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 25 20:43:38 2016
+
+@author: keile
+"""
+import numpy as np
+import unittest
+
+from utils import setmembership
+
+
+class TestUniqueRows(unittest.TestCase):
+    def test_unique_rows_1(self):
+
+        a = np.array([[1, 2], [2, 1], [2, 4], [2, 1], [2, 4]])
+        ua_expected = np.array([[1, 2], [2, 1], [2, 4]])
+        ia_expected = np.array([0, 1, 2])
+        ic_expected = np.array([0, 1, 2, 1, 2])
+        ua, ia, ic = setmembership.unique_rows(a)
+        assert np.sum(np.abs(ua) - np.abs(ua_expected)) == 0
+        assert np.all(ia - ia_expected == 0)
+        assert np.all(ic - ic_expected == 0)
+
+    if __name__ == '__main__':
+        unittest.main()
+
+
+class TestIsmember(unittest.TestCase):
+
+    def test_ismember_rows_with_sort(self):
+        a = np.array([[1, 3, 3, 1, 7],
+                      [3, 3, 2, 3, 0]])
+        b = np.array([[3, 1, 3, 5, 3],
+                      [3, 3, 2, 1, 2]])
+        ma, ia = setmembership.ismember_rows(a, b)
+
+        ma_known = np.array([1, 1, 1, 1, 0], dtype=bool)
+        ia_known = np.array([1, 0, 2, 1])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+
+    def test_ismember_rows_no_sort(self):
+        a = np.array([[1, 3, 3, 1, 7],
+                      [3, 3, 2, 3, 0]])
+        b = np.array([[3, 1, 2, 5, 3],
+                      [3, 3, 3, 1, 1]])
+        ma, ia = setmembership.ismember_rows(a, b, sort=False)
+
+        ma_known = np.array([1, 1, 0, 1, 0], dtype=bool)
+        ia_known = np.array([1, 0, 1])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+    def test_ismember_rows_unqual_sizes_1(self):
+        # a larger than b
+        a = np.array([[1, 3, 3, 1, 7],
+                      [3, 3, 2, 3, 0]])
+        b = np.array([[3, 1, 2, 5],
+                      [3, 3, 3, 1]])
+        ma, ia = setmembership.ismember_rows(a, b)
+
+        ma_known = np.array([1, 1, 1, 1, 0], dtype=bool)
+        ia_known = np.array([1, 0, 2, 1])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+    def test_ismember_rows_unqual_sizes_1(self):
+        # b larger than b
+        a = np.array([[1, 3, 3, 1, 7],
+                      [3, 3, 2, 3, 0]])
+        b = np.array([[3, 1, 2, 5, 3, 4, 7],
+                      [3, 3, 3, 1, 9, 9, 9]])
+        ma, ia = setmembership.ismember_rows(a, b)
+
+        ma_known = np.array([1, 1, 1, 1, 0], dtype=bool)
+        ia_known = np.array([1, 0, 2, 1])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+    def test_ismember_rows_double_occurence_a_no_b(self):
+        # There are duplicate occurences in a that are not found in b
+        a = np.array([[1, 3, 3, 1, 7],
+                      [3, 3, 2, 3, 0]])
+        b = np.array([[3, 2, 5],
+                      [3, 3, 1]])
+        ma, ia = setmembership.ismember_rows(a, b)
+
+        ma_known = np.array([0, 1, 1, 0, 0], dtype=bool)
+        ia_known = np.array([0, 1])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+    def test_ismember_rows_double_occurence_a_and_b(self):
+        # There are duplicate occurences in a, and the same item is found in b
+        a = np.array([[1, 3, 3, 1, 7],
+                      [3, 3, 2, 3, 0]])
+        b = np.array([[3, 1, 2, 5, 3],
+                      [3, 3, 3, 1, 1]])
+        ma, ia = setmembership.ismember_rows(a, b)
+
+        ma_known = np.array([1, 1, 1, 1, 0], dtype=bool)
+        ia_known = np.array([1, 0, 2, 1])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+    def test_ismember_rows_1d(self):
+        a = np.array([0, 2, 1, 3, 0])
+        b = np.array([2, 4, 3])
+
+        ma, ia = setmembership.ismember_rows(a, b)
+
+        ma_known = np.array([0, 1, 0, 1, 0], dtype=bool)
+        ia_known = np.array([0, 2])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+    def test_ismember_rows_1d(self):
+        a = np.array([0, 2, 1, 13, 0])
+        b = np.array([2, 4, 13, 0])
+
+        ma, ia = setmembership.ismember_rows(a, b)
+
+        ma_known = np.array([1, 1, 0, 1, 1], dtype=bool)
+        ia_known = np.array([3, 0, 2, 3])
+
+        assert np.allclose(ma, ma_known)
+        assert np.allclose(ia, ia_known)
+
+    if __name__ == '__main__':
+        unittest.main()
+
+
+class TestUniqueColumns(unittest.TestCase):
+
+    def test_no_common_points(self):
+        p = np.array([[0, 1, 2], [0, 0, 0]])
+        p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p)
+
+        assert np.allclose(p, p_unique)
+        assert np.alltrue(old_2_new == np.arange(3))
+        assert np.alltrue(old_2_new == new_2_old)
+
+    def test_remove_one_point(self):
+        p = np.ones((2, 2))
+        p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p)
+
+        assert np.allclose(p, np.ones((2, 1)))
+        assert np.alltrue(old_2_new == np.zeros(2))
+        assert np.alltrue(new_2_old == np.zeros(1))
+
+    def test_remove_one_of_tree(self):
+        p = np.array([[1, 1, 0], [1, 1, 0]])
+        p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p)
+
+        assert np.allclose(p_unique, np.array([[0, 1], [0, 1]]))
+        assert np.alltrue(old_2_new == np.array([1, 1, 0]))
+        assert np.alltrue(new_2_old == np.array([2, 0]))
+
+    if __name__ == '__main__':
+        unittest.main()
diff --git a/src/porepy/utils/unique.py b/src/porepy/utils/unique.py
new file mode 100644
index 0000000000..0708ce602d
--- /dev/null
+++ b/src/porepy/utils/unique.py
@@ -0,0 +1,124 @@
+import numpy as np
+
+def unique_np1130(ar, return_index=False, return_inverse=False,
+           return_counts=False, axis=None):
+    """
+    Find the unique elements of an array.
+    Returns the sorted unique elements of an array. There are three optional
+    outputs in addition to the unique elements: the indices of the input array
+    that give the unique values, the indices of the unique array that
+    reconstruct the input array, and the number of times each unique value
+    comes up in the input array.
+    Parameters
+    ----------
+    ar : array_like
+        Input array. Unless `axis` is specified, this will be flattened if it
+        is not already 1-D.
+    return_index : bool, optional
+        If True, also return the indices of `ar` (along the specified axis,
+        if provided, or in the flattened array) that result in the unique array.
+    return_inverse : bool, optional
+        If True, also return the indices of the unique array (for the specified
+        axis, if provided) that can be used to reconstruct `ar`.
+    return_counts : bool, optional
+        If True, also return the number of times each unique item appears
+        in `ar`.
+        .. versionadded:: 1.9.0
+    axis : int or None, optional
+        The axis to operate on. If None, `ar` will be flattened beforehand.
+        Otherwise, duplicate items will be removed along the provided axis,
+        with all the other axes belonging to the each of the unique elements.
+        Object arrays or structured arrays that contain objects are not
+        supported if the `axis` kwarg is used.
+        .. versionadded:: 1.13.0
+    Returns
+    -------
+    unique : ndarray
+        The sorted unique values.
+    unique_indices : ndarray, optional
+        The indices of the first occurrences of the unique values in the
+        original array. Only provided if `return_index` is True.
+    unique_inverse : ndarray, optional
+        The indices to reconstruct the original array from the
+        unique array. Only provided if `return_inverse` is True.
+    unique_counts : ndarray, optional
+        The number of times each of the unique values comes up in the
+        original array. Only provided if `return_counts` is True.
+        .. versionadded:: 1.9.0
+    See Also
+    --------
+    numpy.lib.arraysetops : Module with a number of other functions for
+                            performing set operations on arrays.
+    Examples
+    --------
+    >>> np.unique([1, 1, 2, 2, 3, 3])
+    array([1, 2, 3])
+    >>> a = np.array([[1, 1], [2, 3]])
+    >>> np.unique(a)
+    array([1, 2, 3])
+    Return the unique rows of a 2D array
+    >>> a = np.array([[1, 0, 0], [1, 0, 0], [2, 3, 4]])
+    >>> np.unique(a, axis=0)
+    array([[1, 0, 0], [2, 3, 4]])
+    Return the indices of the original array that give the unique values:
+    >>> a = np.array(['a', 'b', 'b', 'c', 'a'])
+    >>> u, indices = np.unique(a, return_index=True)
+    >>> u
+    array(['a', 'b', 'c'],
+           dtype='|S1')
+    >>> indices
+    array([0, 1, 3])
+    >>> a[indices]
+    array(['a', 'b', 'c'],
+           dtype='|S1')
+    Reconstruct the input array from the unique values:
+    >>> a = np.array([1, 2, 6, 4, 2, 3, 2])
+    >>> u, indices = np.unique(a, return_inverse=True)
+    >>> u
+    array([1, 2, 3, 4, 6])
+    >>> indices
+    array([0, 1, 4, 3, 1, 2, 1])
+    >>> u[indices]
+    array([1, 2, 6, 4, 2, 3, 2])
+    """
+    ar = np.asanyarray(ar)
+    if axis is None:
+        return np.unique(ar, return_index, return_inverse, return_counts)
+    if not (-ar.ndim <= axis < ar.ndim):
+        raise ValueError('Invalid axis kwarg specified for unique')
+
+    ar = np.swapaxes(ar, axis, 0)
+    orig_shape, orig_dtype = ar.shape, ar.dtype
+    # Must reshape to a contiguous 2D array for this to work...
+    ar = ar.reshape(orig_shape[0], -1)
+    ar = np.ascontiguousarray(ar)
+
+    if ar.dtype.char in (np.typecodes['AllInteger'] +
+                         np.typecodes['Datetime'] + 'S'):
+        # Optimization: Creating a view of your data with a np.void data type of
+        # size the number of bytes in a full row. Handles any type where items
+        # have a unique binary representation, i.e. 0 is only 0, not +0 and -0.
+        dtype = np.dtype((np.void, ar.dtype.itemsize * ar.shape[1]))
+    else:
+        dtype = [('f{i}'.format(i=i), ar.dtype) for i in range(ar.shape[1])]
+
+    try:
+        consolidated = ar.view(dtype)
+    except TypeError:
+        # There's no good way to do this for object arrays, etc...
+        msg = 'The axis argument to unique is not supported for dtype {dt}'
+        raise TypeError(msg.format(dt=ar.dtype))
+
+    def reshape_uniq(uniq):
+        uniq = uniq.view(orig_dtype)
+        uniq = uniq.reshape(-1, *orig_shape[1:])
+        uniq = np.swapaxes(uniq, 0, axis)
+        return uniq
+
+    output = np.unique(consolidated, return_index,
+                       return_inverse, return_counts)
+    if not (return_index or return_inverse or return_counts):
+        return reshape_uniq(output)
+    else:
+        uniq = reshape_uniq(output[0])
+        return (uniq,) + output[1:]