Add utils to new common repository

.gitignore

@@ -2,3 +2,4 @@
 *~
 *.py~
 *.swp
+*.py#

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

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
+
+#------------------------------------------------------------------------------#

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

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]
+
+#------------------------------------------------------------------------------#