Python FieldVariable.create_output方法代码示例

# 需要导入模块: from sfepy.discrete import FieldVariable [as 别名]
# 或者: from sfepy.discrete.FieldVariable import create_output [as 别名]
def save_basis_on_mesh(mesh, options, output_dir, lin,
                       permutations=None, suffix=''):
    if permutations is not None:
        mesh = mesh.copy()
        gel = GeometryElement(mesh.descs[0])
        perms = gel.get_conn_permutations()[permutations]
        conn = mesh.cmesh.get_cell_conn()
        n_el, n_ep = conn.num, gel.n_vertex
        offsets = nm.arange(n_el) * n_ep

        conn.indices[:] = conn.indices.take((perms + offsets[:, None]).ravel())

    domain = FEDomain('domain', mesh)

    omega = domain.create_region('Omega', 'all')
    field = Field.from_args('f', nm.float64, shape=1, region=omega,
                            approx_order=options.max_order,
                            poly_space_base=options.basis)
    var = FieldVariable('u', 'unknown', field)

    if options.plot_dofs:
        import sfepy.postprocess.plot_dofs as pd
        import sfepy.postprocess.plot_cmesh as pc
        ax = pc.plot_wireframe(None, mesh.cmesh)
        ax = pd.plot_global_dofs(ax, field.get_coor(), field.econn)
        ax = pd.plot_local_dofs(ax, field.get_coor(), field.econn)
        if options.dofs is not None:
            ax = pd.plot_nodes(ax, field.get_coor(), field.econn,
                               field.poly_space.nodes,
                               get_dofs(options.dofs, var.n_dof))
        pd.plt.show()

    output('dofs: %d' % var.n_dof)

    vec = nm.empty(var.n_dof, dtype=var.dtype)
    n_digit, _format = get_print_info(var.n_dof, fill='0')
    name_template = os.path.join(output_dir,
                                 'dof_%s%s.vtk' % (_format, suffix))
    for ip in get_dofs(options.dofs, var.n_dof):
        output('dof %d...' % ip)

        vec.fill(0.0)
        vec[ip] = 1.0

        var.set_data(vec)

        if options.derivative == 0:
            out = var.create_output(vec, linearization=lin)

        else:
            out = create_expression_output('ev_grad.ie.Elements(u)',
                                           'u', 'f', {'f' : field}, None,
                                           Variables([var]),
                                           mode='qp', verbose=False,
                                           min_level=lin.min_level,
                                           max_level=lin.max_level,
                                           eps=lin.eps)

        name = name_template % ip
        ensure_path(name)
        out['u'].mesh.write(name, out=out)

        output('...done (%s)' % name)

# 需要导入模块: from sfepy.discrete import FieldVariable [as 别名]
# 或者: from sfepy.discrete.FieldVariable import create_output [as 别名]

#.........这里部分代码省略.........
                                              verbose=True)

    output('...done in', time.clock() - tt)

    output('creating solver...')
    tt = time.clock()

    conf = Struct(method='bcgsl', precond='jacobi', sub_precond=None,
                  i_max=10000, eps_a=1e-50, eps_r=1e-6, eps_d=1e4,
                  verbose=True)
    status = {}
    ls = PETScKrylovSolver(conf, comm=comm, mtx=pmtx, status=status)

    field_ranges = {}
    for ii, variable in enumerate(variables.iter_state(ordered=True)):
        field_ranges[variable.name] = lfds[ii].petsc_dofs_range

    ls.set_field_split(field_ranges, comm=comm)

    ev = PETScParallelEvaluator(pb, pdofs, drange, True,
                                psol, comm, verbose=True)

    nls_status = {}
    conf = Struct(method='newtonls',
                  i_max=5, eps_a=0, eps_r=1e-5, eps_s=0.0,
                  verbose=True)
    nls = PETScNonlinearSolver(conf, pmtx=pmtx, prhs=prhs, comm=comm,
                               fun=ev.eval_residual,
                               fun_grad=ev.eval_tangent_matrix,
                               lin_solver=ls, status=nls_status)

    output('...done in', time.clock() - tt)

    output('solving...')
    tt = time.clock()

    state = pb.create_state()
    state.apply_ebc()

    ev.psol_i[...] = state()
    ev.gather(psol, ev.psol_i)

    psol = nls(psol)

    ev.scatter(ev.psol_i, psol)
    sol0_i = ev.psol_i[...]

    output('...done in', time.clock() - tt)

    output('saving solution...')
    tt = time.clock()

    state.set_full(sol0_i)
    out = state.create_output_dict()

    filename = os.path.join(options.output_dir, 'sol_%02d.h5' % comm.rank)
    pb.domain.mesh.write(filename, io='auto', out=out)

    gather_to_zero = pl.create_gather_to_zero(psol)

    psol_full = gather_to_zero(psol)

    if comm.rank == 0:
        sol = psol_full[...].copy()

        u = FieldVariable('u', 'parameter', field1,
                          primary_var_name='(set-to-None)')
        remap = gfds[0].id_map
        ug = sol[remap]

        p = FieldVariable('p', 'parameter', field2,
                          primary_var_name='(set-to-None)')
        remap = gfds[1].id_map
        pg = sol[remap]

        if (((order_u == 1) and (order_p == 1))
            or (options.linearization == 'strip')):
            out = u.create_output(ug)
            out.update(p.create_output(pg))
            filename = os.path.join(options.output_dir, 'sol.h5')
            mesh.write(filename, io='auto', out=out)

        else:
            out = u.create_output(ug, linearization=Struct(kind='adaptive',
                                                           min_level=0,
                                                           max_level=order_u,
                                                           eps=1e-3))

            filename = os.path.join(options.output_dir, 'sol_u.h5')
            out['u'].mesh.write(filename, io='auto', out=out)

            out = p.create_output(pg, linearization=Struct(kind='adaptive',
                                                           min_level=0,
                                                           max_level=order_p,
                                                           eps=1e-3))

            filename = os.path.join(options.output_dir, 'sol_p.h5')
            out['p'].mesh.write(filename, io='auto', out=out)

    output('...done in', time.clock() - tt)

# 需要导入模块: from sfepy.discrete import FieldVariable [as 别名]
# 或者: from sfepy.discrete.FieldVariable import create_output [as 别名]
    def test_continuity(self):
        from sfepy.base.base import Struct
        from sfepy.mesh.mesh_generators import gen_block_mesh
        from sfepy.discrete import FieldVariable
        from sfepy.discrete.fem import FEDomain, Field
        from sfepy.discrete.projections import make_l2_projection_data
        import sfepy.discrete.fem.refine_hanging as rh

        dims = [1.5, 2.0, 1.3]
        shape = [3, 3, 3]
        centre = [0.0, 0.0, 0.0]
        probe_gens = {"2_4": _gen_lines_2_4, "3_8": _gen_grid_3_8}

        ok = True
        for key in self.gel_names:
            gel = self.gels[key]
            probe_gen = probe_gens[key]

            perms = gel.get_conn_permutations()

            dim = gel.dim
            for io, order in enumerate(range(1, 4)):
                mesh00 = gen_block_mesh(dims[:dim], shape[:dim], centre[:dim], name="block")

                for ip, perm in enumerate(perms):
                    self.report("geometry: %s, order: %d, permutation: %d: %s" % (key, order, ip, perm))
                    mesh0 = mesh00.copy()
                    conn = mesh0.cmesh.get_conn(dim, 0).indices
                    conn = conn.reshape((mesh0.n_el, -1))
                    conn[-1, :] = conn[-1, perm]

                    domain0 = FEDomain("d", mesh0)

                    refine = nm.zeros(mesh0.n_el, dtype=nm.uint8)
                    refine[:-1] = 1

                    subs = None
                    domain, subs = rh.refine(domain0, refine, subs=subs)

                    omega = domain.create_region("Omega", "all")
                    field = Field.from_args("fu", nm.float64, 1, omega, approx_order=order)

                    field.substitute_dofs(subs)

                    uvar = FieldVariable("u", "parameter", field, primary_var_name="(set-to-None)")

                    field.restore_dofs(store=True)
                    field.substitute_dofs(subs=None, restore=True)

                    make_l2_projection_data(uvar, eval_fun)

                    field.restore_dofs()

                    bbox = domain.get_mesh_bounding_box()
                    eps = 1e-7

                    save = False
                    for ii, (probe0, probe1) in enumerate(probe_gen(bbox, eps)):
                        probe0.set_options(close_limit=0.0)
                        probe1.set_options(close_limit=0.0)

                        pars0, vals0 = probe0(uvar)
                        pars1, vals1 = probe1(uvar)

                        assert_(nm.allclose(pars0, pars1, atol=1e-14, rtol=0.0))

                        _ok = nm.allclose(vals0, vals1, atol=20.0 * eps, rtol=0.0)
                        if not _ok:
                            save = True
                            self.report("probe %d failed! (max. error: %e)" % (ii, nm.abs(vals0 - vals1).max()))

                        ok = ok and _ok

                    if (ip == 0) or save:
                        out = uvar.create_output()
                        filenames = _build_filenames(self.options.out_dir, key, order, ip)

                        domain.mesh.write(filenames[0], out=out)

                        linearization = Struct(kind="adaptive", min_level=0, max_level=4, eps=1e-2)

                        out = uvar.create_output(linearization=linearization)
                        val = out["u"]
                        mesh = val.get("mesh", domain.mesh)
                        mesh.write(filenames[1], out=out)

        return ok