Python Domain.setup_facets方法代码示例

# 需要导入模块: from sfepy.fem import Domain [as 别名]
# 或者: from sfepy.fem.Domain import setup_facets [as 别名]
def main():
    parser = OptionParser(usage = usage, version = "%prog " + sfepy.__version__)
    parser.add_option( "-m", "--mesh",
                       action = "store_true", dest = "save_mesh",
                       default = True,
                       help = "save surface mesh [default: %default]" )
    parser.add_option( "-n", "--no-surface",
                       action = "store_true", dest = "no_surface",
                       default = False,
                       help = "do not output surface [default: %default]" )
    (options, args) = parser.parse_args()

    if (len( args ) == 2):
        filename_in = args[0];
        filename_out = args[1];
    else:
        parser.print_help(),
        return

    if (filename_in == '-'):
        file_in = sys.stdin
    else:
        file_in = open( filename_in, "r" ); 

    mesh = Mesh.from_file( filename_in )

    if (filename_in != '-'):
        file_in.close()

    domain = Domain('domain', mesh)
    domain.setup_groups()

    if domain.has_faces():
        domain.fix_element_orientation()
        domain.setup_facets(create_edges=False)

        lst, surf_faces = domain.surface_faces()

        surf_mesh = Mesh.from_surface( surf_faces, mesh )

        if options.save_mesh:
            aux = edit_filename(filename_in, prefix='surf_', new_ext='.mesh')
            surf_mesh.write(aux, io='auto')

        if options.no_surface:
            return

        gr_s = surface_graph( surf_faces, mesh.n_nod )
##         import sfepy.base.plotutils as plu
##         plu.spy( gr_s )
##         plu.pylab.show()

        n_comp, comps = surface_components( gr_s, surf_faces )
#        print 'components:', n_comp

        ccs, comps = comps, nm.zeros( (0,1), nm.int32 )
        for cc in ccs:
            comps = nm.concatenate( (comps, cc[:,nm.newaxis]), 0 )

        out = nm.concatenate( (lst, comps), 1 )

        if (filename_out == '-'):
            file_out = sys.stdout
        else:
            file_out = open( filename_out, "w" ); 
        for row in out:
            file_out.write( '%d %d %d %d\n' % (row[0], row[1], row[2], row[3]) )
        if (filename_out != '-'):
            file_out.close()

# 需要导入模块: from sfepy.fem import Domain [as 别名]
# 或者: from sfepy.fem.Domain import setup_facets [as 别名]

#.........这里部分代码省略.........
        for ip in get_dofs(options.dofs, ps.n_nod):
            output('shape function %d...' % ip)

            def eval_dofs(iels, rx):
                if options.derivative == 0:
                    bf = ps.eval_base(rx).squeeze()
                    rvals = bf[None, :, ip:ip+1]

                else:
                    bfg = ps.eval_base(rx, diff=True)
                    rvals = bfg[None, ..., ip]

                return rvals

            def eval_coors(iels, rx):
                bf = gps.eval_base(rx).squeeze()
                coors = nm.dot(bf, gel.coors)[None, ...]
                return coors

            (level, coors, conn,
             vdofs, mat_ids) = create_output(eval_dofs, eval_coors, 1,
                                             ps, min_level=lin.min_level,
                                             max_level=lin.max_level,
                                             eps=lin.eps)
            out = {
                'bf' : Struct(name='output_data',
                              mode='vertex', data=vdofs,
                              var_name='bf', dofs=None)
            }

            mesh = Mesh.from_data('bf_mesh', coors, None, [conn], [mat_ids],
                                  [options.geometry])

            name = name_template % ip
            mesh.write(name, out=out)

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

    else:
        mesh = Mesh.from_file(options.mesh)
        output('mesh geometry:')
        output('  dimension: %d, vertices: %d, elements: %d'
               % (mesh.dim, mesh.n_nod, mesh.n_el))

        domain = Domain('domain', mesh)

        if options.permutations:
            permutations = [int(ii) for ii in options.permutations.split(',')]
            output('using connectivity permutations:', permutations)
            for group in domain.iter_groups():
                perms = group.gel.get_conn_permutations()[permutations]
                offsets = nm.arange(group.shape.n_el) * group.shape.n_ep

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

            domain.setup_facets()

        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, 1)

        if options.plot_dofs:
            import sfepy.postprocess.plot_dofs as pd
            group = domain.groups[0]
            ax = pd.plot_mesh(None, mesh.coors, mesh.conns[0], group.gel.edges)
            ax = pd.plot_global_dofs(ax, field.get_coor(), field.aps[0].econn)
            ax = pd.plot_local_dofs(ax, field.get_coor(), field.aps[0].econn)
            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.vtk' % _format)
        for ip in get_dofs(options.dofs, var.n_dof):
            output('dof %d...' % ip)

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

            var.data_from_any(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
            out['u'].mesh.write(name, out=out)

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

# 需要导入模块: from sfepy.fem import Domain [as 别名]
# 或者: from sfepy.fem.Domain import setup_facets [as 别名]
def save_basis_on_mesh(mesh, options, output_dir, lin,
                       permutations=None, suffix=''):
    domain = Domain('domain', mesh)

    if permutations is not None:
        for group in domain.iter_groups():
            perms = group.gel.get_conn_permutations()[permutations]
            offsets = nm.arange(group.shape.n_el) * group.shape.n_ep

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

        domain.setup_facets()

    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, 1)

    if options.plot_dofs:
        import sfepy.postprocess.plot_dofs as pd
        group = domain.groups[0]
        ax = pd.plot_mesh(None, mesh.coors, mesh.conns[0], group.gel.edges)
        ax = pd.plot_global_dofs(ax, field.get_coor(), field.aps[0].econn)
        ax = pd.plot_local_dofs(ax, field.get_coor(), field.aps[0].econn)
        if options.dofs is not None:
            ax = pd.plot_nodes(ax, field.get_coor(), field.aps[0].econn,
                               field.aps[0].interp.poly_spaces['v'].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)