Python Mesh.from_region方法代码示例

# 需要导入模块: from sfepy.discrete.fem import Mesh [as 别名]
# 或者: from sfepy.discrete.fem.Mesh import from_region [as 别名]
def main():
    parser = ArgumentParser(description=__doc__.rstrip(),
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('filename', help=helps['filename'])
    parser.add_argument('-d', '--detailed',
                        action='store_true', dest='detailed',
                        default=False, help=helps['detailed'])
    options = parser.parse_args()

    mesh = Mesh.from_file(options.filename)

    output(mesh.cmesh)
    output('element types:', mesh.descs)
    output('nodal BCs:', sorted(mesh.nodal_bcs.keys()))

    bbox = mesh.get_bounding_box()
    output('bounding box:\n%s'
           % '\n'.join('%s: [%14.7e, %14.7e]' % (name, bbox[0, ii], bbox[1, ii])
                       for ii, name in enumerate('xyz'[:mesh.dim])))

    output('centre:           [%s]'
           % ', '.join('%14.7e' % ii for ii in 0.5 * (bbox[0] + bbox[1])))
    output('coordinates mean: [%s]'
           % ', '.join('%14.7e' % ii for ii in mesh.coors.mean(0)))

    if not options.detailed: return

    domain = FEDomain(mesh.name, mesh)

    for dim in range(1, mesh.cmesh.tdim + 1):
        volumes = mesh.cmesh.get_volumes(dim)
        output('volumes of %d %dD entities:\nmin: %.7e mean: %.7e median:'
               ' %.7e max: %.7e'
               % (mesh.cmesh.num[dim], dim, volumes.min(), volumes.mean(),
                  nm.median(volumes), volumes.max()))

    euler = lambda mesh: nm.dot(mesh.cmesh.num, [1, -1, 1, -1])
    ec = euler(mesh)
    output('Euler characteristic:', ec)

    graph = mesh.create_conn_graph(verbose=False)
    n_comp, _ = graph_components(graph.shape[0], graph.indptr, graph.indices)
    output('number of connected components:', n_comp)

    if mesh.dim > 1:
        region = domain.create_region('surf', 'vertices of surface', 'facet')
        surf_mesh = Mesh.from_region(region, mesh,
                                     localize=True, is_surface=True)
        FEDomain(surf_mesh.name, surf_mesh) # Calls CMesh.setup_entities().

        sec = euler(surf_mesh)
        output('surface Euler characteristic:', sec)
        if mesh.dim == 3:
            output('surface genus:', (2.0 - sec) / 2.0)

        surf_graph = surf_mesh.create_conn_graph(verbose=False)
        n_comp, _ = graph_components(surf_graph.shape[0],
                                     surf_graph.indptr, surf_graph.indices)
        output('number of connected surface components:', n_comp)

# 需要导入模块: from sfepy.discrete.fem import Mesh [as 别名]
# 或者: from sfepy.discrete.fem.Mesh import from_region [as 别名]
def refine_region(domain0, region0, region1):
    """
    Coarse cell sub_cells[ii, 0] in mesh0 is split into sub_cells[ii, 1:] in
    mesh1.

    The new fine cells are interleaved among the original coarse cells so that
    the indices of the coarse cells do not change.

    The cell groups are preserved. The vertex groups are preserved only in the
    coarse (non-refined) cells.
    """
    if region1 is None:
        return domain0, None

    mesh0 = domain0.mesh
    mesh1 = Mesh.from_region(region1, mesh0)
    domain1 = FEDomain('d', mesh1)
    domain1r = domain1.refine()
    mesh1r = domain1r.mesh

    n_cell = region1.shape.n_cell
    n_sub = 4 if mesh0.cmesh.tdim == 2 else 8

    sub_cells = nm.empty((n_cell, n_sub + 1), dtype=nm.uint32)
    sub_cells[:, 0] = region1.cells
    sub_cells[:, 1] = region1.cells
    aux = nm.arange((n_sub - 1) * n_cell, dtype=nm.uint32)
    sub_cells[:, 2:] = mesh0.n_el + aux.reshape((n_cell, -1))

    coors0, vgs0, conns0, mat_ids0, descs0 = mesh0._get_io_data()
    coors, vgs, _conns, _mat_ids, descs = mesh1r._get_io_data()

    # Preserve vertex groups of non-refined cells.
    vgs[:len(vgs0)] = vgs0

    def _interleave_refined(c0, c1):
        if c1.ndim == 1:
            c0 = c0[:, None]
            c1 = c1[:, None]

        n_row, n_col = c1.shape
        n_new = region0.shape.n_cell + n_row

        out = nm.empty((n_new, n_col), dtype=c0.dtype)
        out[region0.cells] = c0[region0.cells]
        out[region1.cells] = c1[::n_sub]
        aux = c1.reshape((-1, n_col * n_sub))
        out[mesh0.n_el:] = aux[:, n_col:].reshape((-1, n_col))

        return out

    conn = _interleave_refined(conns0[0], _conns[0])
    mat_id = _interleave_refined(mat_ids0[0], _mat_ids[0]).squeeze()

    mesh = Mesh.from_data('a', coors, vgs, [conn], [mat_id], descs)
    domain = FEDomain('d', mesh)

    return domain, sub_cells

# 需要导入模块: from sfepy.discrete.fem import Mesh [as 别名]
# 或者: from sfepy.discrete.fem.Mesh import from_region [as 别名]
def create_local_problem(omega_gi, orders):
    """
    Local problem definition using a domain corresponding to the global region
    `omega_gi`.
    """
    order_u, order_p = orders

    mesh = omega_gi.domain.mesh

    # All tasks have the whole mesh.
    bbox = mesh.get_bounding_box()
    min_x, max_x = bbox[:, 0]
    eps_x = 1e-8 * (max_x - min_x)

    min_y, max_y = bbox[:, 1]
    eps_y = 1e-8 * (max_y - min_y)

    mesh_i = Mesh.from_region(omega_gi, mesh, localize=True)
    domain_i = FEDomain('domain_i', mesh_i)
    omega_i = domain_i.create_region('Omega', 'all')

    gamma1_i = domain_i.create_region('Gamma1',
                                      'vertices in (x < %.10f)'
                                      % (min_x + eps_x),
                                      'facet', allow_empty=True)
    gamma2_i = domain_i.create_region('Gamma2',
                                      'vertices in (x > %.10f)'
                                      % (max_x - eps_x),
                                      'facet', allow_empty=True)
    gamma3_i = domain_i.create_region('Gamma3',
                                      'vertices in (y < %.10f)'
                                      % (min_y + eps_y),
                                      'facet', allow_empty=True)

    field1_i = Field.from_args('fu', nm.float64, mesh.dim, omega_i,
                               approx_order=order_u)

    field2_i = Field.from_args('fp', nm.float64, 1, omega_i,
                               approx_order=order_p)

    output('field 1: number of local DOFs:', field1_i.n_nod)
    output('field 2: number of local DOFs:', field2_i.n_nod)

    u_i = FieldVariable('u_i', 'unknown', field1_i, order=0)
    v_i = FieldVariable('v_i', 'test', field1_i, primary_var_name='u_i')
    p_i = FieldVariable('p_i', 'unknown', field2_i, order=1)
    q_i = FieldVariable('q_i', 'test', field2_i, primary_var_name='p_i')

    if mesh.dim == 2:
        alpha = 1e2 * nm.array([[0.132], [0.132], [0.092]])

    else:
        alpha = 1e2 * nm.array([[0.132], [0.132], [0.132],
                                [0.092], [0.092], [0.092]])

    mat = Material('m', D=stiffness_from_lame(mesh.dim, lam=10, mu=5),
                   k=1, alpha=alpha)
    integral = Integral('i', order=2*(max(order_u, order_p)))

    t11 = Term.new('dw_lin_elastic(m.D, v_i, u_i)',
                   integral, omega_i, m=mat, v_i=v_i, u_i=u_i)
    t12 = Term.new('dw_biot(m.alpha, v_i, p_i)',
                   integral, omega_i, m=mat, v_i=v_i, p_i=p_i)
    t21 = Term.new('dw_biot(m.alpha, u_i, q_i)',
                   integral, omega_i, m=mat, u_i=u_i, q_i=q_i)
    t22 = Term.new('dw_laplace(m.k, q_i, p_i)',
                   integral, omega_i, m=mat, q_i=q_i, p_i=p_i)

    eq1 = Equation('eq1', t11 - t12)
    eq2 = Equation('eq1', t21 + t22)
    eqs = Equations([eq1, eq2])

    ebc1 = EssentialBC('ebc1', gamma1_i, {'u_i.all' : 0.0})
    ebc2 = EssentialBC('ebc2', gamma2_i, {'u_i.0' : 0.05})
    def bc_fun(ts, coors, **kwargs):
        val = 0.3 * nm.sin(4 * nm.pi * (coors[:, 0] - min_x) / (max_x - min_x))
        return val

    fun = Function('bc_fun', bc_fun)
    ebc3 = EssentialBC('ebc3', gamma3_i, {'p_i.all' : fun})

    pb = Problem('problem_i', equations=eqs, active_only=False)
    pb.time_update(ebcs=Conditions([ebc1, ebc2, ebc3]))
    pb.update_materials()

    return pb

# 需要导入模块: from sfepy.discrete.fem import Mesh [as 别名]
# 或者: from sfepy.discrete.fem.Mesh import from_region [as 别名]
def main():
    parser = ArgumentParser(description=__doc__,
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('-s', '--scale', metavar='scale',
                        action='store', dest='scale',
                        default=None, help=helps['scale'])
    parser.add_argument('-c', '--center', metavar='center',
                        action='store', dest='center',
                        default=None, help=helps['center'])
    parser.add_argument('-r', '--refine', metavar='level',
                        action='store', type=int, dest='refine',
                        default=0, help=helps['refine'])
    parser.add_argument('-f', '--format', metavar='format',
                        action='store', type=str, dest='format',
                        default=None, help=helps['format'])
    parser.add_argument('-l', '--list', action='store_true',
                        dest='list', help=helps['list'])
    parser.add_argument('-m', '--merge', action='store_true',
                        dest='merge', help=helps['merge'])
    parser.add_argument('-t', '--tri-tetra', action='store_true',
                        dest='tri_tetra', help=helps['tri-tetra'])
    parser.add_argument('-2', '--2d', action='store_true',
                        dest='force_2d', help=helps['2d'])
    parser.add_argument('--save-per-mat', action='store_true',
                        dest='save_per_mat', help=helps['save-per-mat'])
    parser.add_argument('--remesh', metavar='options',
                        action='store', dest='remesh',
                        default=None, help=helps['remesh'])
    parser.add_argument('filename_in')
    parser.add_argument('filename_out')
    options = parser.parse_args()

    if options.list:
        output('Supported readable mesh formats:')
        output('--------------------------------')
        output_mesh_formats('r')
        output('')
        output('Supported writable mesh formats:')
        output('--------------------------------')
        output_mesh_formats('w')
        sys.exit(0)

    scale = _parse_val_or_vec(options.scale, 'scale', parser)
    center = _parse_val_or_vec(options.center, 'center', parser)

    filename_in = options.filename_in
    filename_out = options.filename_out

    if options.remesh:
        import tempfile
        import shlex
        import subprocess

        dirname = tempfile.mkdtemp()

        is_surface = options.remesh.startswith('q')
        if is_surface:
            mesh = Mesh.from_file(filename_in)
            domain = FEDomain(mesh.name, mesh)
            region = domain.create_region('surf', 'vertices of surface',
                                          'facet')
            surf_mesh = Mesh.from_region(region, mesh,
                                         localize=True, is_surface=True)

            filename = op.join(dirname, 'surf.mesh')
            surf_mesh.write(filename, io='auto')

        else:
            import shutil

            shutil.copy(filename_in, dirname)
            filename = op.join(dirname, op.basename(filename_in))

        qopts = ''.join(options.remesh.split()) # Remove spaces.
        command = 'tetgen -BFENkACp%s %s' % (qopts, filename)
        args = shlex.split(command)
        subprocess.call(args)

        root, ext = op.splitext(filename)
        mesh = Mesh.from_file(root + '.1.vtk')

        remove_files(dirname)

    else:
        mesh = Mesh.from_file(filename_in)

    if options.force_2d:
        data = list(mesh._get_io_data())
        data[0] = data[0][:, :2]
        mesh = Mesh.from_data(mesh.name, *data)

    if scale is not None:
        if len(scale) == 1:
            tr = nm.eye(mesh.dim, dtype=nm.float64) * scale
        elif len(scale) == mesh.dim:
            tr = nm.diag(scale)
        else:
            raise ValueError('bad scale! (%s)' % scale)
        mesh.transform_coors(tr)

#.........这里部分代码省略.........

# 需要导入模块: from sfepy.discrete.fem import Mesh [as 别名]
# 或者: from sfepy.discrete.fem.Mesh import from_region [as 别名]
def main():
    parser = OptionParser(usage=usage, version="%prog " + sfepy.__version__)
    parser.add_option("-m", "--mesh",
                      action="store_true", dest="save_mesh",
                      default=False,
                      help="save surface mesh")
    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 = FEDomain('domain', mesh)

    if options.save_mesh:
        region = domain.create_region('surf', 'vertices of surface', 'facet')
        surf_mesh = Mesh.from_region(region, mesh,
                                     localize=True, is_surface=True)
        aux = edit_filename(filename_in, prefix='surf_', new_ext='.mesh')
        surf_mesh.write(aux, io='auto')

    if domain.has_faces():
        domain.fix_element_orientation()

        lst, surf_faces = get_surface_faces(domain)

        if options.no_surface:
            return

        gr_s = surface_graph(surf_faces, mesh.n_nod)

        n_comp, comps = surface_components(gr_s, surf_faces)
        output('number of surface 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\n' % (row[0], row[1], row[2]))
        if (filename_out != '-'):
            file_out.close()

# 需要导入模块: from sfepy.discrete.fem import Mesh [as 别名]
# 或者: from sfepy.discrete.fem.Mesh import from_region [as 别名]
def main():
    parser = ArgumentParser(description=__doc__,
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument("--version", action="version",
                        version="%(prog)s " + sfepy.__version__)
    parser.add_argument("-m", "--mesh",
                        action="store_true", dest="save_mesh",
                        default=False,
                        help="save surface mesh")
    parser.add_argument("-n", "--no-surface",
                        action="store_true", dest="no_surface",
                        default=False,
                        help="do not output surface [default: %(default)s]")
    parser.add_argument('filename_in', help="'-' is for stdin")
    parser.add_argument('filename_out', help="'-' is for  stdout")
    options = parser.parse_args()

    filename_in = options.filename_in
    filename_out = options.filename_out

    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 = FEDomain('domain', mesh)

    if options.save_mesh:
        region = domain.create_region('surf', 'vertices of surface', 'facet')
        surf_mesh = Mesh.from_region(region, mesh,
                                     localize=True, is_surface=True)
        aux = edit_filename(filename_in, prefix='surf_', new_ext='.mesh')
        surf_mesh.write(aux, io='auto')

    if domain.has_faces():
        domain.fix_element_orientation()

        lst, surf_faces = get_surface_faces(domain)

        if options.no_surface:
            return

        gr_s = surface_graph(surf_faces, mesh.n_nod)

        n_comp, comps = surface_components(gr_s, surf_faces)
        output('number of surface 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\n' % (row[0], row[1], row[2]))
        if (filename_out != '-'):
            file_out.close()

# 需要导入模块: from sfepy.discrete.fem import Mesh [as 别名]
# 或者: from sfepy.discrete.fem.Mesh import from_region [as 别名]
def create_local_problem(omega_gi, orders):
    """
    Local problem definition using a domain corresponding to the global region
    `omega_gi`.
    """
    order_u, order_p = orders

    mesh = omega_gi.domain.mesh

    # All tasks have the whole mesh.
    bbox = mesh.get_bounding_box()
    min_x, max_x = bbox[:, 0]
    eps_x = 1e-8 * (max_x - min_x)

    min_y, max_y = bbox[:, 1]
    eps_y = 1e-8 * (max_y - min_y)

    mesh_i = Mesh.from_region(omega_gi, mesh, localize=True)
    domain_i = FEDomain("domain_i", mesh_i)
    omega_i = domain_i.create_region("Omega", "all")

    gamma1_i = domain_i.create_region("Gamma1", "vertices in (x < %.10f)" % (min_x + eps_x), "facet", allow_empty=True)
    gamma2_i = domain_i.create_region("Gamma2", "vertices in (x > %.10f)" % (max_x - eps_x), "facet", allow_empty=True)
    gamma3_i = domain_i.create_region("Gamma3", "vertices in (y < %.10f)" % (min_y + eps_y), "facet", allow_empty=True)

    field1_i = Field.from_args("fu", nm.float64, mesh.dim, omega_i, approx_order=order_u)

    field2_i = Field.from_args("fp", nm.float64, 1, omega_i, approx_order=order_p)

    output("field 1: number of local DOFs:", field1_i.n_nod)
    output("field 2: number of local DOFs:", field2_i.n_nod)

    u_i = FieldVariable("u_i", "unknown", field1_i, order=0)
    v_i = FieldVariable("v_i", "test", field1_i, primary_var_name="u_i")
    p_i = FieldVariable("p_i", "unknown", field2_i, order=1)
    q_i = FieldVariable("q_i", "test", field2_i, primary_var_name="p_i")

    if mesh.dim == 2:
        alpha = 1e2 * nm.array([[0.132], [0.132], [0.092]])

    else:
        alpha = 1e2 * nm.array([[0.132], [0.132], [0.132], [0.092], [0.092], [0.092]])

    mat = Material("m", lam=10, mu=5, k=1, alpha=alpha)
    integral = Integral("i", order=2 * (max(order_u, order_p)))

    t11 = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v_i, u_i)", integral, omega_i, m=mat, v_i=v_i, u_i=u_i)
    t12 = Term.new("dw_biot(m.alpha, v_i, p_i)", integral, omega_i, m=mat, v_i=v_i, p_i=p_i)
    t21 = Term.new("dw_biot(m.alpha, u_i, q_i)", integral, omega_i, m=mat, u_i=u_i, q_i=q_i)
    t22 = Term.new("dw_laplace(m.k, q_i, p_i)", integral, omega_i, m=mat, q_i=q_i, p_i=p_i)

    eq1 = Equation("eq1", t11 - t12)
    eq2 = Equation("eq1", t21 + t22)
    eqs = Equations([eq1, eq2])

    ebc1 = EssentialBC("ebc1", gamma1_i, {"u_i.all": 0.0})
    ebc2 = EssentialBC("ebc2", gamma2_i, {"u_i.0": 0.05})

    def bc_fun(ts, coors, **kwargs):
        val = 0.3 * nm.sin(4 * nm.pi * (coors[:, 0] - min_x) / (max_x - min_x))
        return val

    fun = Function("bc_fun", bc_fun)
    ebc3 = EssentialBC("ebc3", gamma3_i, {"p_i.all": fun})

    pb = Problem("problem_i", equations=eqs, active_only=False)
    pb.time_update(ebcs=Conditions([ebc1, ebc2, ebc3]))
    pb.update_materials()

    return pb