def get_edge_graph(self):
        """
        Return the graph of region edges as a sparse matrix having uid(k) + 1
        at (i, j) if vertex[i] is connected with vertex[j] by the edge k.

        Degenerate edges are ignored.
        """
        from scipy.sparse import csr_matrix

        ed = self.domain.ed

        rows, cols, vals = [], [], []
        for ig, edges in self.edges.iteritems():
            e_verts = ed.facets[edges]
            ii = nm.where(e_verts[:, 0] != e_verts[:, 1])[0]
            edges = edges[ii]
            e_verts = e_verts[ii]

            vals.append(ed.uid_i[edges] + 1)
            rows.append(e_verts[:, 0])
            cols.append(e_verts[:, 1])

        vals, indx = nm.unique(nm.concatenate(vals), return_index=True)
        rows = nm.concatenate(rows)[indx]
        cols = nm.concatenate(cols)[indx]

        num = self.all_vertices.max() + 1
        graph = csr_matrix((vals, (rows, cols)), shape=(num, num))

        nnz = graph.nnz
        # Symmetrize.
        graph = graph + graph.T
        assert_(graph.nnz == 2 * nnz)

        return graph

def norm_l2_along_axis(ar, axis=1, n_item=None, squared=False):
    """Compute l2 norm of rows (axis=1) or columns (axis=0) of a 2D array.

    n_item ... use only the first n_item columns/rows
    squared ... if True, return the norm squared
    """
    assert_(axis in [0, 1])
    assert_(ar.ndim == 2)

    other = 1 - axis
    vec = nm.zeros((ar.shape[other],), dtype=nm.float64)

    if n_item is None:
        n_item = ar.shape[axis]
    else:
        n_item = min( n_item, ar.shape[axis] )

    if axis == 1:
        for ii in range( n_item ):
            vec += ar[:,ii]**2
    else:
        for ii in range( n_item ):
            vec += ar[ii,:]**2

    if not squared:
        vec = nm.sqrt( vec )

    return vec

    def set_dofs(self, fun=0.0, region=None, dpn=None, warn=None):
        """
        Set the values of DOFs in a given region using a function of space
        coordinates or value `fun`.
        """
        if region is None:
            region = self.region

        if dpn is None:
            dpn = self.shape[0]

        aux = self.get_dofs_in_region(region, clean=True, warn=warn)
        nods = nm.unique(nm.hstack(aux))

        if callable(fun):
            vals = fun(self.get_coor(nods))

        elif nm.isscalar(fun):
            vals = nm.repeat([fun], nods.shape[0] * dpn)

        elif isinstance(fun, nm.ndarray):
            assert_(len(fun) == dpn)
            vals = nm.repeat(fun, nods.shape[0])

        else:
            raise ValueError('unknown function/value type! (%s)' % type(fun))

        return nods, vals

def barycentric_coors(coors, s_coors):
    """
    Get barycentric (area in 2D, volume in 3D) coordinates of points
    with coordinates `coors` w.r.t. the simplex given by `s_coors`.

    Returns
    -------
    bc : array
        The barycentric coordinates. Then reference element coordinates
        `xi = dot(bc.T, ref_coors)`.
    """
    n_v, dim = s_coors.shape
    n_c, dim2 = coors.shape
    assert_(dim == dim2)
    assert_(n_v == (dim + 1))

    mtx = nm.ones((n_v, n_v), nm.float64)
    mtx[0:n_v-1,:] = s_coors.T

    rhs = nm.empty((n_v,n_c), nm.float64)
    rhs[0:n_v-1,:] = coors.T
    rhs[n_v-1,:] = 1.0

    bc = nla.solve(mtx, rhs)

    return bc

    def get_edge_graph(self):
        """
        Return the graph of region edges as a sparse matrix having uid(k) + 1
        at (i, j) if vertex[i] is connected with vertex[j] by the edge k.

        Degenerate edges are ignored.
        """
        from scipy.sparse import csr_matrix

        cmesh = self.domain.cmesh

        e_verts = cmesh.get_incident(0, self.edges, 1)
        e_verts.shape = (e_verts.shape[0] / 2, 2)

        ii = nm.where(e_verts[:, 0] != e_verts[:, 1])[0]
        edges = self.edges[ii]
        e_verts = e_verts[ii]

        vals = edges + 1
        rows = e_verts[:, 0]
        cols = e_verts[:, 1]

        num = self.vertices.max() + 1
        graph = csr_matrix((vals, (rows, cols)), shape=(num, num))

        nnz = graph.nnz
        # Symmetrize.
        graph = graph + graph.T
        assert_(graph.nnz == 2 * nnz)

        return graph

    def __init__(self, anchor, normal, bounds):
        Struct.__init__(self, anchor=nm.array(anchor, dtype=nm.float64),
                        bounds=nm.asarray(bounds, dtype=nm.float64))
        self.normal = nm.asarray(normal, dtype=nm.float64)

        norm = nm.linalg.norm
        self.normal /= norm(self.normal)

        e3 = [0.0, 0.0, 1.0]
        dd = nm.dot(e3, self.normal)
        rot_angle = nm.arccos(dd)

        if nm.abs(rot_angle) < 1e-14:
            mtx = nm.eye(3, dtype=nm.float64)
            bounds2d = self.bounds[:, :2]

        else:
            rot_axis = nm.cross([0.0, 0.0, 1.0], self.normal)
            mtx = la.make_axis_rotation_matrix(rot_axis, rot_angle)

            mm = la.insert_strided_axis(mtx, 0, self.bounds.shape[0])
            rbounds = la.dot_sequences(mm, self.bounds)
            bounds2d = rbounds[:, :2]

        assert_(nm.allclose(nm.dot(mtx, self.normal), e3,
                            rtol=0.0, atol=1e-12))

        self.adotn = nm.dot(self.anchor, self.normal)

        self.rot_angle = rot_angle
        self.mtx = mtx
        self.bounds2d = bounds2d

    def __call__(self, problem=None, data=None):
        problem = get_default(problem, self.problem)

        problem.set_equations(self.equations)

        problem.select_bcs(ebc_names=self.ebcs, epbc_names=self.epbcs,
                           lcbc_names=self.get_default_attr('lcbcs', []))

        problem.update_materials(problem.ts)

        self.init_solvers(problem)

        variables = problem.get_variables()

        if hasattr(self, 'set_variables'):
            if isinstance(self.set_variables, list):
                self.set_variables_default(variables, self.set_variables,
                                           data)
            else:
                self.set_variables(variables, **data)

        state = problem.solve()
        assert_(state.has_ebc())

        corr_sol = CorrSolution(name=self.name,
                                state=state.get_parts())

        self.save(corr_sol, problem)

        return corr_sol

    def _petsc_call(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                    i_max=None, mtx=None, status=None, comm=None,
                    context=None, **kwargs):
        tt = time.clock()

        conf = get_default(conf, self.conf)
        mtx = get_default(mtx, self.mtx)
        status = get_default(status, self.status)
        context = get_default(context, self.context)
        comm = get_default(comm, self.comm)

        mshape = mtx.size if isinstance(mtx, self.petsc.Mat) else mtx.shape
        rshape = [rhs.size] if isinstance(rhs, self.petsc.Vec) else rhs.shape

        assert_(mshape[0] == mshape[1] == rshape[0])
        if x0 is not None:
            xshape = [x0.size] if isinstance(x0, self.petsc.Vec) else x0.shape
            assert_(xshape[0] == rshape[0])

        result = call(self, rhs, x0, conf, eps_a, eps_r, i_max, mtx, status,
                      comm, context=context, **kwargs)

        ttt = time.clock() - tt
        if status is not None:
            status['time'] = ttt
            status['n_iter'] = self.ksp.getIterationNumber()

        return result

    def set_dofs(self, fun=0.0, region=None, dpn=None, warn=None):
        """
        Set the values of given DOFs using a function of space coordinates or
        value `fun`.

        Notes
        -----
        Works for a constant value over an entire patch side only.
        """
        if region is None:
            region = self.region

        if dpn is None:
            dpn = self.n_components

        nods = []
        vals = []

        aux = self.get_dofs_in_region(region, clean=True, warn=warn)
        nods = nm.unique(nm.hstack(aux))

        if nm.isscalar(fun):
            vals = nm.repeat([fun], nods.shape[0] * dpn)

        elif isinstance(fun, nm.ndarray):
            assert_(len(fun) == dpn)
            vals = nm.repeat(fun, nods.shape[0])

        else:
            raise ValueError('unknown function/value type! (%s)' % type(fun))

        return nods, vals

def insert_sparse_to_csr(mtx1, mtx2, irs, ics):
    """
    Insert a sparse matrix `mtx2` into a CSR sparse matrix `mtx1` at
    rows `irs` and columns `ics`. The submatrix `mtx1[irs,ics]` must
    already be preallocated and have the same structure as `mtx2`.
    """
    import sfepy.discrete.common.extmods.assemble as asm

    if isinstance(irs, slice):
        irs = nm.arange(irs.start, irs.stop, irs.step, dtype=nm.int32)

    if isinstance(ics, slice):
        ics = nm.arange(ics.start, ics.stop, ics.step, dtype=nm.int32)

    n_row, n_col = mtx1.shape

    assert_((irs.min() >= 0) and (irs.max() < n_row))
    assert_((ics.min() >= 0) and (ics.max() < n_col))

    aux = mtx2.tocoo()
    data = nm.ascontiguousarray(aux.data[:,None,None,None])
    rows = irs[aux.row[:,None]]
    cols = ics[aux.col[:,None]]

    iels = nm.arange(rows.shape[0], dtype=nm.int32)
    asm.assemble_matrix(mtx1.data, mtx1.indptr, mtx1.indices, data,
                        iels, 1.0, rows, cols)

def compute_bezier_extraction(knots, degrees):
    """
    Compute local (element) Bezier extraction operators for a nD B-spline
    parametric domain.

    Parameters
    ----------
    knots : sequence of array or array
        The knot vectors.
    degrees : sequence of ints or int
        Polynomial degrees in each parametric dimension.

    Returns
    -------
    cs : list of lists of 2D arrays
        The element extraction operators in each parametric dimension.
    """
    if isinstance(degrees, int): degrees = [degrees]

    knots = _get_knots_tuple(knots)
    dim = len(knots)
    assert_(dim == len(degrees))

    cs = []
    for ii, knots1d in enumerate(knots):
        cs1d = compute_bezier_extraction_1d(knots1d, degrees[ii])
        cs.append(cs1d)

    return cs

def get_poly(order, dim, is_simplex=False):
    """
    Construct a polynomial of given `order` in space dimension `dim`,
    and integrate it symbolically over a rectangular or simplex domain
    for coordinates in [0, 1].
    """
    xs = symarray("x", dim)

    opd = max(1, int((order + 1) / dim))

    poly = 1.0
    oo = 0
    for ii, x in enumerate(xs):
        if ((oo + opd) > order) or (ii == (len(xs) - 1)):
            opd = max(order - oo, 0)

        poly *= x ** opd + 1
        oo += opd

    assert_(oo == order)

    limits = [[xs[ii], 0, 1] for ii in range(dim)]
    if is_simplex:
        for ii in range(1, dim):
            for ip in range(0, ii):
                limits[ii][2] -= xs[ip]

    integral = sm.integrate(poly, *reversed(limits))

    return xs, poly, limits, integral

 def __call__(self, parser, namespace, value, option_string=None):
     vals = value.split(',')
     assert_(len(vals) in [2, 3, 6])
     val = tuple(float(ii) for ii in vals)
     if len(vals) == 6:
         val = val[:3] + (list(val[3:]),)
     setattr(namespace, self.dest, val)

    def __call__( self, problem = None, data = None, save_hook = None ):
        """data: corrs_pressure, evp, optionally vec_g"""
        problem = get_default( problem, self.problem )
        ts = problem.get_time_solver().ts

        corrs, evp = [data[ii] for ii in self.requires[:2]]
        if len(self.requires) == 3:
            vec_g = data[self.requires[2]]
        else:
            vec_g = None

        assert_( evp.ebcs == self.ebcs )
        assert_( evp.epbcs == self.epbcs )

        filename = self.get_dump_name()
        savename = self.get_save_name()

        self.setup_equations(self.equations)

        solve = self.compute_correctors
        solve(evp, 1.0, corrs.state, ts, filename, savename, vec_g=vec_g)

        if self.check:
            self.setup_equations(self.verify_equations)
            self.init_solvers(problem)

            output( 'verifying correctors %s...' % self.name )
            verify = self.verify_correctors
            ok = verify(1.0, corrs.state, filename)
            output( '...done, ok: %s' % ok )

        return Struct( name = self.name,
                       filename = filename )

    def _standard_call(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                       i_max=None, mtx=None, status=None, context=None,
                       **kwargs):
        tt = time.clock()

        conf = get_default(conf, self.conf)
        mtx = get_default(mtx, self.mtx)
        status = get_default(status, self.status)
        context = get_default(context, self.context)

        assert_(mtx.shape[0] == mtx.shape[1] == rhs.shape[0])
        if x0 is not None:
            assert_(x0.shape[0] == rhs.shape[0])

        result = call(self, rhs, x0, conf, eps_a, eps_r, i_max, mtx, status,
                      context=context, **kwargs)
        if isinstance(result, tuple):
            result, n_iter = result

        else:
            n_iter = -1 # Number of iterations is undefined/unavailable.

        ttt = time.clock() - tt
        if status is not None:
            status['time'] = ttt
            status['n_iter'] = n_iter

        return result

    def __call__(self, volume=None, problem=None, data=None):
        problem = get_default(problem, self.problem)
        opts = self.app_options

        iw_dir = nm.array(opts.incident_wave_dir, dtype=nm.float64)
        dim = problem.get_dim()
        assert_(dim == iw_dir.shape[0])

        iw_dir = iw_dir / nla.norm(iw_dir)

        dispersion = data[self.requires[0]]

        wave_vectors = dispersion.logs.eig_vectors

        pas = []

        iw_dir = iw_dir / nla.norm(iw_dir)
        idims = range(iw_dir.shape[0])
        pi2 = 0.5 * nm.pi
        for vecs in wave_vectors:
            pa = nm.empty(vecs.shape[:-1], dtype=nm.float64)
            for ir, vec in enumerate(vecs):
                for ic in idims:
                    vv = vec[:,ic]
                    # Ensure the angle is in [0, pi/2].
                    val = nm.arccos(nm.dot(iw_dir, vv) / nla.norm(vv))
                    if val > pi2:
                        val = nm.pi - val
                    pa[ir,ic] = val

            pas.append(pa)

        return pas

        def get_constants(ts=None, coors=None, mode=None,
                          term=None, problem=None, **kwargs):
            out = {}
            if mode == 'qp':
                qps = term.get_physical_qps()
                assert_(qps.num == coors.shape[0])

                for key, val in six.iteritems(values):
                    if '.' in key: continue
                    rval = nm.array(val[list(val.keys())[0]], dtype=nm.float64,
                                    ndmin=3)
                    s0 = rval.shape[1:]
                    matdata = nm.zeros(qps.shape[:2] + s0, dtype=nm.float64)

                    for rkey, rval in six.iteritems(val):
                        region = problem.domain.regions[rkey]
                        rval = nm.array(rval, dtype=nm.float64, ndmin=3)

                        cells = region.get_cells(true_cells_only=False)
                        ii = term.region.get_cell_indices(cells,
                                                          true_cells_only=False)
                        matdata[ii] = rval

                    out[key] = matdata.reshape((-1,) + s0)

            return out

def postprocess(filename_input, filename_results, options):
    """
    Postprocess probe data files - replot, integrate data.
    """
    from matplotlib import pyplot as plt

    header, results = read_results(filename_input,
                                   only_names=options.only_names)
    output(header)

    fig = plt.figure()
    for name, result in results.iteritems():
        pars, vals = result[:, 0], result[:, 1]

        ii = nm.where(nm.isfinite(vals))[0]
        # Nans only at the edges.
        assert_(nm.diff(ii).sum() == (len(ii)-1))

        val = integrate_along_line(pars[ii], vals[ii], options.radial)

        label = r'%s: $\int\ %s' % (name, name)
        if options.radial:
            label += ' (r)'
        label += '$ = %.5e'% val

        plt.plot(pars, vals, label=label, lw=0.2, marker='+', ms=1)
        plt.ylabel('probed data')
        plt.xlabel('probe coordinate')

        output(label)

    plt.legend()

    fig.savefig(filename_results)

    def solve_pressure_eigenproblem(self, mtx, eig_problem=None,
                                    n_eigs=0, check=False):
        """G = B*AI*BT or B*AI*BT+D"""

        def get_slice(n_eigs, nn):
            if n_eigs > 0:
                ii = slice(0, n_eigs)
            elif n_eigs < 0:
                ii = slice(nn + n_eigs, nn)
            else:
                ii = slice(0, 0)
            return ii

        eig_problem = get_default(eig_problem, self.eig_problem)
        n_eigs = get_default(n_eigs, self.n_eigs)
        check = get_default(check, self.check)

        mtx_c, mtx_b, action_aibt = mtx['C'], mtx['B'], mtx['action_aibt']
        mtx_g = mtx_b * action_aibt.to_array() # mtx_b must be sparse!
        if eig_problem == 'B*AI*BT+D':
            mtx_g += mtx['D'].toarray()

        mtx['G'] = mtx_g
        output(mtx_c.shape, mtx_g.shape)

        eigs, mtx_q = eig(mtx_c.toarray(), mtx_g, method='eig.sgscipy')

        if check:
            ee = nm.diag(sc.dot(mtx_q.T * mtx_c, mtx_q)).squeeze()
            oo = nm.diag(sc.dot(sc.dot(mtx_q.T,  mtx_g), mtx_q)).squeeze()
            try:
                assert_(nm.allclose(ee, eigs))
                assert_(nm.allclose(oo, nm.ones_like(eigs)))
            except ValueError:
                debug()

        nn = mtx_c.shape[0]
        if isinstance(n_eigs, tuple):
            output('required number of eigenvalues: (%d, %d)' % n_eigs)
            if sum(n_eigs) < nn:
                ii0 = get_slice(n_eigs[0], nn)
                ii1 = get_slice(-n_eigs[1], nn)
                eigs = nm.concatenate((eigs[ii0], eigs[ii1]))
                mtx_q = nm.concatenate((mtx_q[:,ii0], mtx_q[:,ii1]), 1) 
        else:
            output('required number of eigenvalues: %d' % n_eigs)
            if (n_eigs != 0) and (abs(n_eigs) < nn):
                ii = get_slice(n_eigs, nn)
                eigs = eigs[ii]
                mtx_q = mtx_q[:,ii]

##         from sfepy.base.plotutils import pylab, iplot
##         pylab.semilogy(eigs)
##         pylab.figure(2)
##         iplot(eigs)
##         pylab.show()
##         debug()

        out = Struct(eigs=eigs, mtx_q=mtx_q)
        return out

    def _setup_vertex_dofs(self):
        """
        Setup vertex DOF connectivity.
        """
        if self.node_desc.vertex is None:
            return 0, None

        region = self.region

        cmesh = self.domain.cmesh
        conn, offsets = cmesh.get_incident(0, region.cells, region.tdim,
                                           ret_offsets=True)

        vertices = nm.unique(conn)
        remap = prepare_remap(vertices, region.n_v_max)
        n_dof = vertices.shape[0]

        aux = nm.unique(nm.diff(offsets))
        assert_(len(aux) == 1, 'region with multiple reference geometries!')
        offset = aux[0]


        # Remap vertex node connectivity to field-local numbering.
        aux = conn.reshape((-1, offset)).astype(nm.int32)
        self.econn[:, :offset] = nm.take(remap, aux)

        return n_dof, remap