Python base.assert_函数代码示例

    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