Python Mesh.from_data方法代码示例

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
    def save_axes(self, filename):
        coors = []
        conns = []
        mat_ids = []
        offset = 0
        for ig, dual_surface in self.dual_surfaces.iteritems():
            cc = nm.r_[dual_surface.edge_centre_coors,
                       dual_surface.dual_coors]
            coors.append(cc)

            conn = dual_surface.conn.copy() + offset
            conn[:,1:] += dual_surface.edge_centre_coors.shape[0]
            conns.append(conn)
            
            mat_id = nm.empty((conn.shape[0],), dtype=nm.int32)
            mat_id[:] = ig
            mat_ids.append(mat_id)

            offset += cc.shape[0]

        coors = nm.concatenate(coors, axis=0)

        out = {}
        for ig, dual_surface in self.dual_surfaces.iteritems():
            eto = edge_data_to_output
            out['en_%d' % ig] = eto(coors, conns[ig], dual_surface.e_sort,
                                    dual_surface.edge_normals)
            out['ed_%d' % ig] = eto(coors, conns[ig], dual_surface.e_sort,
                                    dual_surface.edge_dirs)
            out['eo_%d' % ig] = eto(coors, conns[ig], dual_surface.e_sort,
                                    dual_surface.edge_ortho)

        dual_mesh = Mesh.from_data('dual_mesh_vectors', coors, None, conns,
                                   mat_ids, ['2_4'] * len(conns))
        dual_mesh.write(filename, io='auto', out=out)

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def refine_3_8(mesh_in, ed, fa):
    """
    Refines hexahedral mesh by cutting cutting each edge in half and
    making 8 new finer hexahedrons out of one coarser one.
    """
    # Unique edge centres.
    e_coors, e_uid = ed.get_coors()
    e_centres = 0.5 * nm.sum(e_coors, axis=1)

    # Unique face centres.
    f_coors, f_uid = fa.get_coors()
    f_centres = 0.25 * nm.sum(f_coors, axis=1)

    # Unique element centres.
    coors = mesh_in.get_element_coors()
    centres = 0.125 * nm.sum(coors, axis=1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres, f_centres, centres]

    o1 = mesh_in.n_nod
    o2 = o1 + e_centres.shape[0]
    o3 = o2 + f_centres.shape[0]

    st = nm.vstack

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        e_indx = ed.indx[ig]
        f_indx = fa.indx[ig]
        off = mesh_in.el_offsets[ig]
        n_el  = conn.shape[0]

        e_nodes = ed.uid_i[e_indx].reshape((n_el, 12)) + o1
        f_nodes = fa.uid_i[f_indx].reshape((n_el, 6)) + o2
        nodes = nm.arange(n_el) + off + o3

        c = nm.c_[conn, e_nodes, f_nodes, nodes].T

        new_conn = st([c[0], c[8], c[20], c[11], c[16], c[22], c[26], c[21],
                       c[1], c[9], c[20], c[8], c[17], c[24], c[26], c[22],
                       c[2], c[10], c[20], c[9], c[18], c[25], c[26], c[24],
                       c[3], c[11], c[20], c[10], c[19], c[21], c[26], c[25],
                       c[4], c[15], c[23], c[12], c[16], c[21], c[26], c[22],
                       c[5], c[12], c[23], c[13], c[17], c[22], c[26], c[24],
                       c[6], c[13], c[23], c[14], c[18], c[24], c[26], c[25],
                       c[7], c[14], c[23], c[15], c[19], c[25], c[26], c[21]]).T
        new_conn = new_conn.reshape((8 * n_el, 8))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(8)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, conns,
                          mat_ids, mesh_in.descs )

    return mesh

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def refine_3_4(mesh_in, ed):
    """
    Refines tetrahedra by cutting each edge in half and making 8 new
    finer tetrahedra out of one coarser one. Old nodal coordinates come
    first in `coors`, then the new ones. The new tetrahedra are similar
    to the old one, no degeneration is supposed to occur as at most 3
    congruence classes of tetrahedra appear, even when re-applied
    iteratively (provided that `conns` are not modified between two
    applications - ordering of vertices in tetrahedra matters not only
    for positivity of volumes).

    References:

    - Juergen Bey: Simplicial grid refinement: on Freudenthal s algorithm and 
      the optimal number of congruence classes, Numer.Math. 85 (2000), 
      no. 1, 1--29, or
    - Juergen Bey: Tetrahedral grid refinement, Computing 55 (1995), 
      no. 4, 355--378, or
      http://citeseer.ist.psu.edu/bey95tetrahedral.html
    """
    # Unique edge centres.
    e_coors, e_uid = ed.get_coors()
    e_centres = 0.5 * nm.sum(e_coors, axis=1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres]

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        indx = ed.indx[ig]
        n_el  = conn.shape[0]

        e_nodes = ed.uid_i[indx].reshape((n_el, 6)) + mesh_in.n_nod

        c = nm.c_[conn, e_nodes].T

        new_conn = nm.vstack([c[0], c[4], c[6], c[7],
                              c[4], c[1], c[5], c[8],
                              c[6], c[5], c[2], c[9],
                              c[7], c[8], c[9], c[3],
                              c[4], c[6], c[7], c[8],
                              c[4], c[6], c[8], c[5],
                              c[6], c[7], c[8], c[9],
                              c[6], c[5], c[9], c[8]]).T
        new_conn = new_conn.reshape((8 * n_el, 4))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(8)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, conns,
                          mat_ids, mesh_in.descs )

    return mesh

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def refine_2_4(mesh_in, ed):
    """
    Refines mesh out of quadrilaterals by cutting cutting each edge in
    half and making 4 new finer quadrilaterals out of one coarser one.
    """
    # Unique edge centres.
    e_coors, e_uid = ed.get_coors()
    e_centres = 0.5 * nm.sum(e_coors, axis=1)

    # Unique element centres.
    coors = mesh_in.get_element_coors()
    centres = 0.25 * nm.sum(coors, axis=1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres, centres]

    o1 = mesh_in.n_nod
    o2 = o1 + e_centres.shape[0]

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        e_indx = ed.indx[ig]
        off = mesh_in.el_offsets[ig]
        n_el  = conn.shape[0]

        e_nodes = ed.uid_i[e_indx].reshape((n_el, 4)) + o1
        nodes = nm.arange(n_el) + off + o2

        c = nm.c_[conn, e_nodes, nodes].T

        new_conn = nm.vstack([c[0], c[4], c[8], c[7],
                              c[1], c[5], c[8], c[4],
                              c[2], c[6], c[8], c[5],
                              c[3], c[7], c[8], c[6]]).T
        new_conn = new_conn.reshape((4 * n_el, 4))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(4)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, conns,
                          mat_ids, mesh_in.descs )

    return mesh

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def refine_2_4(mesh_in, cmesh):
    """
    Refines mesh out of quadrilaterals by cutting cutting each edge in
    half and making 4 new finer quadrilaterals out of one coarser one.
    """
    # Unique edge centres.
    e_centres = cmesh.get_centroids(cmesh.dim - 1)

    # Unique element centres.
    centres = cmesh.get_centroids(cmesh.dim)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres, centres]

    o1 = mesh_in.n_nod
    o2 = o1 + e_centres.shape[0]

    cc = cmesh.get_conn(cmesh.dim, cmesh.dim - 1)
    offs = cc.offsets

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        off0, off1 = mesh_in.el_offsets[ig : ig + 2]
        n_el = conn.shape[0]

        e_nodes = cc.indices[offs[off0] : offs[off1]].reshape((n_el, 4)) + o1
        nodes = nm.arange(n_el) + off0 + o2

        c = nm.c_[conn, e_nodes, nodes].T

        new_conn = nm.vstack(
            [c[0], c[4], c[8], c[7], c[1], c[5], c[8], c[4], c[2], c[6], c[8], c[5], c[3], c[7], c[8], c[6]]
        ).T
        new_conn = new_conn.reshape((4 * n_el, 4))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(4)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + "_r", coors, None, conns, mat_ids, mesh_in.descs)

    return mesh

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
    def save(self, filename):
        coors = []
        conns = []
        mat_ids = []
        offset = 0
        for ig, dual_surface in self.dual_surfaces.iteritems():
            cc = dual_surface.dual_coors
            coors.append(cc)

            conn = dual_surface.conn[:, 1:].copy() + offset
            conns.append(conn)

            mat_id = nm.empty((conn.shape[0],), dtype=nm.int32)
            mat_id[:] = ig
            mat_ids.append(mat_id)

            offset += cc.shape[0]

        coors = nm.concatenate(coors, axis=0)

        dual_mesh = Mesh.from_data("dual_mesh", coors, None, conns, mat_ids, ["2_3"] * len(conns))
        dual_mesh.write(filename, io="auto")

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def refine_2_3(mesh_in, ed):
    """
    Refines mesh out of triangles by cutting cutting each edge in half
    and making 4 new finer triangles out of one coarser one.
    """
    # Unique edge centres.
    e_coors, e_uid = ed.get_coors()
    e_centres = 0.5 * nm.sum(e_coors, axis=1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres]

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        indx = ed.indx[ig]
        n_el  = conn.shape[0]

        e_nodes = ed.uid_i[indx].reshape((n_el, 3)) + mesh_in.n_nod

        c = nm.c_[conn, e_nodes].T

        new_conn = nm.vstack([c[0], c[3], c[5],
                              c[3], c[4], c[5],
                              c[1], c[4], c[3],
                              c[2], c[5], c[4]]).T
        new_conn = new_conn.reshape((4 * n_el, 3))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(4)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, conns,
                          mat_ids, mesh_in.descs )

    return mesh

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

#.........这里部分代码省略.........
    etype : integer, optional
        'q' - quadrilateral or hexahedral elements
        't' - triangular or tetrahedral elements
    Returns
    -------
    mesh : mesh
        Finite element mesh.
    """

    dims = dims.squeeze()
    dim = len(dims)
    nddims = nm.array(voxels.shape) + 1

    nodemtx = nm.zeros(nddims, dtype=nm.int32)

    if dim == 2:
        #iy, ix = nm.where(voxels.transpose())
        iy, ix = nm.where(voxels)
        nel = ix.shape[0]

        if etype == 'q':
            nodemtx[ix,iy] += 1
            nodemtx[ix + 1,iy] += 1
            nodemtx[ix + 1,iy + 1] += 1
            nodemtx[ix,iy + 1] += 1

        elif etype == 't':
            nodemtx[ix,iy] += 2
            nodemtx[ix + 1,iy] += 1
            nodemtx[ix + 1,iy + 1] += 2
            nodemtx[ix,iy + 1] += 1
            nel *= 2

    elif dim == 3:
        #iy, ix, iz = nm.where(voxels.transpose(1, 0, 2))
        iy, ix, iz = nm.where(voxels)
        nel = ix.shape[0]

        if etype == 'q':
            nodemtx[ix,iy,iz] += 1
            nodemtx[ix + 1,iy,iz] += 1
            nodemtx[ix + 1,iy + 1,iz] += 1
            nodemtx[ix,iy + 1,iz] += 1
            nodemtx[ix,iy,iz + 1] += 1
            nodemtx[ix + 1,iy,iz + 1] += 1
            nodemtx[ix + 1,iy + 1,iz + 1] += 1
            nodemtx[ix,iy + 1,iz + 1] += 1

        elif etype == 't':
            nodemtx[ix,iy,iz] += 6
            nodemtx[ix + 1,iy,iz] += 2
            nodemtx[ix + 1,iy + 1,iz] += 2
            nodemtx[ix,iy + 1,iz] += 2
            nodemtx[ix,iy,iz + 1] += 2
            nodemtx[ix + 1,iy,iz + 1] += 2
            nodemtx[ix + 1,iy + 1,iz + 1] += 6
            nodemtx[ix,iy + 1,iz + 1] += 2
            nel *= 6

    else:
        msg = 'incorrect voxel dimension! (%d)' % dim
        raise ValueError(msg)

    ndidx = nm.where(nodemtx)
    coors = nm.array(ndidx).transpose() * dims
    nnod = coors.shape[0]

    nodeid = -nm.ones(nddims, dtype=nm.int32)
    nodeid[ndidx] = nm.arange(nnod)

    # generate elements
    if dim == 2:
        elems = nm.array([nodeid[ix,iy],
                          nodeid[ix + 1,iy],
                          nodeid[ix + 1,iy + 1],
                          nodeid[ix,iy + 1]]).transpose()

    elif dim == 3:
        elems = nm.array([nodeid[ix,iy,iz],
                          nodeid[ix + 1,iy,iz],
                          nodeid[ix + 1,iy + 1,iz],
                          nodeid[ix,iy + 1,iz],
                          nodeid[ix,iy,iz + 1],
                          nodeid[ix + 1,iy,iz + 1],
                          nodeid[ix + 1,iy + 1,iz + 1],
                          nodeid[ix,iy + 1,iz + 1]]).transpose()

    if etype == 't':
        elems = elems_q2t(elems)

    eid = etype + str(dim)
    eltab = {'q2': 4, 'q3': 8, 't2': 3, 't3': 4}

    mesh = Mesh.from_data('voxel_data',
                          coors, nm.ones((nnod, ), dtype=nm.int32),
                          {0: nm.ascontiguousarray(elems)},
                          {0: nm.ones((nel, ), dtype=nm.int32)},
                          {0: '%d_%d' % (dim, eltab[eid])})

    return mesh

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def main():
    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-b', '--basis', metavar='name',
                      action='store', dest='basis',
                      default='lagrange', help=help['basis'])
    parser.add_option('-d', '--derivative', metavar='d', type=int,
                      action='store', dest='derivative',
                      default=0, help=help['derivative'])
    parser.add_option('-n', '--max-order', metavar='order', type=int,
                      action='store', dest='max_order',
                      default=2, help=help['max_order'])
    parser.add_option('-g', '--geometry', metavar='name',
                      action='store', dest='geometry',
                      default='2_4', help=help['geometry'])
    parser.add_option('-m', '--mesh', metavar='mesh',
                      action='store', dest='mesh',
                      default=None, help=help['mesh'])
    parser.add_option('', '--permutations', metavar='permutations',
                      action='store', dest='permutations',
                      default=None, help=help['permutations'])
    parser.add_option('', '--dofs', metavar='dofs',
                      action='store', dest='dofs',
                      default=None, help=help['dofs'])
    parser.add_option('-l', '--lin-options', metavar='options',
                      action='store', dest='lin_options',
                      default='min_level=2,max_level=5,eps=1e-3',
                      help=help['lin_options'])
    parser.add_option('', '--plot-dofs',
                      action='store_true', dest='plot_dofs',
                      default=False, help=help['plot_dofs'])
    options, args = parser.parse_args()

    if len(args) == 1:
        output_dir = args[0]
    else:
        parser.print_help(),
        return

    output('polynomial space:', options.basis)
    output('max. order:', options.max_order)

    lin = Struct(kind='adaptive', min_level=2, max_level=5, eps=1e-3)
    for opt in options.lin_options.split(','):
        key, val = opt.split('=')
        setattr(lin, key, eval(val))

    if options.mesh is None:
        dim, n_ep = int(options.geometry[0]), int(options.geometry[2])
        output('reference element geometry:')
        output('  dimension: %d, vertices: %d' % (dim, n_ep))

        gel = GeometryElement(options.geometry)
        gps = PolySpace.any_from_args(None, gel, 1,
                                      base=options.basis)
        ps = PolySpace.any_from_args(None, gel, options.max_order,
                                     base=options.basis)

        n_digit, _format = get_print_info(ps.n_nod, fill='0')
        name_template = os.path.join(output_dir, 'bf_%s.vtk' % _format)
        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:')
#.........这里部分代码省略.........

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def main():
    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-b', '--basis', metavar='name',
                      action='store', dest='basis',
                      default='lagrange', help=help['basis'])
    parser.add_option('-d', '--derivative', metavar='d', type=int,
                      action='store', dest='derivative',
                      default=0, help=help['derivative'])
    parser.add_option('-n', '--max-order', metavar='order', type=int,
                      action='store', dest='max_order',
                      default=2, help=help['max_order'])
    parser.add_option('-g', '--geometry', metavar='name',
                      action='store', dest='geometry',
                      default='2_4', help=help['geometry'])
    parser.add_option('-m', '--mesh', metavar='mesh',
                      action='store', dest='mesh',
                      default=None, help=help['mesh'])
    parser.add_option('', '--permutations', metavar='permutations',
                      action='store', dest='permutations',
                      default=None, help=help['permutations'])
    parser.add_option('', '--dofs', metavar='dofs',
                      action='store', dest='dofs',
                      default=None, help=help['dofs'])
    parser.add_option('-l', '--lin-options', metavar='options',
                      action='store', dest='lin_options',
                      default='min_level=2,max_level=5,eps=1e-3',
                      help=help['lin_options'])
    parser.add_option('', '--plot-dofs',
                      action='store_true', dest='plot_dofs',
                      default=False, help=help['plot_dofs'])
    options, args = parser.parse_args()

    if len(args) == 1:
        output_dir = args[0]
    else:
        parser.print_help(),
        return

    output('polynomial space:', options.basis)
    output('max. order:', options.max_order)

    lin = Struct(kind='adaptive', min_level=2, max_level=5, eps=1e-3)
    for opt in options.lin_options.split(','):
        key, val = opt.split('=')
        setattr(lin, key, eval(val))

    if options.mesh is None:
        dim, n_ep = int(options.geometry[0]), int(options.geometry[2])
        output('reference element geometry:')
        output('  dimension: %d, vertices: %d' % (dim, n_ep))

        gel = GeometryElement(options.geometry)
        gps = PolySpace.any_from_args(None, gel, 1,
                                      base=options.basis)
        ps = PolySpace.any_from_args(None, gel, options.max_order,
                                     base=options.basis)

        n_digit, _format = get_print_info(ps.n_nod, fill='0')
        name_template = os.path.join(output_dir, 'bf_%s.vtk' % _format)
        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
            ensure_path(name)
            mesh.write(name, out=out)

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

    else:
        mesh = Mesh.from_file(options.mesh)
#.........这里部分代码省略.........

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def main():
    parser = OptionParser(usage=usage, version="%prog")
    parser.add_option(
        "-o", "", metavar="filename", action="store", dest="output_filename", default="out.vtk", help=help["filename"]
    )
    parser.add_option(
        "-d", "--dims", metavar="dims", action="store", dest="dims", default="[1.0, 1.0, 1.0]", help=help["dims"]
    )
    parser.add_option(
        "-s", "--shape", metavar="shape", action="store", dest="shape", default="[11, 11, 11]", help=help["shape"]
    )
    parser.add_option(
        "-c",
        "--centre",
        metavar="centre",
        action="store",
        dest="centre",
        default="[0.0, 0.0, 0.0]",
        help=help["centre"],
    )
    (options, args) = parser.parse_args()

    dims = eval("nm.array( %s, dtype = nm.float64 )" % options.dims)
    shape = eval("nm.array( %s, dtype = nm.int32 )" % options.shape)
    centre = eval("nm.array( %s, dtype = nm.float64 )" % options.centre)

    print dims
    print shape
    print centre

    dim = shape.shape[0]

    x0 = centre - 0.5 * dims
    dd = dims / (shape - 1)

    grid = nm.zeros(shape, dtype=nm.float64)
    n_nod = nm.prod(shape)
    coors = nm.zeros((n_nod, dim + 1), dtype=nm.float64)

    # This is 3D only...
    bar = MyBar("       nodes:")
    bar.init(n_nod)
    for ii, ic in enumerate(cycle(shape)):
        ix, iy, iz = ic
        grid[ix, iy, iz] = ii
        coors[ii, :-1] = x0 + ic * dd
        if not (ii % 100):
            bar.update(ii)
    print
    n_el = nm.prod(shape - 1)
    conn = nm.zeros((n_el, 8), dtype=nm.int32)
    bar = MyBar("       elements:")
    bar.init(n_el)
    for ii, (ix, iy, iz) in enumerate(cycle(shape - 1)):
        conn[ii, :] = [
            grid[ix, iy, iz],
            grid[ix + 1, iy, iz],
            grid[ix + 1, iy + 1, iz],
            grid[ix, iy + 1, iz],
            grid[ix, iy, iz + 1],
            grid[ix + 1, iy, iz + 1],
            grid[ix + 1, iy + 1, iz + 1],
            grid[ix, iy + 1, iz + 1],
        ]
        if not (ii % 100):
            bar.update(ii)
    print
    mat_id = nm.zeros((n_el,), dtype=nm.int32)
    desc = "3_8"

    mesh = Mesh.from_data(options.output_filename, coors, [conn], [mat_id], [desc])
    mesh.write(options.output_filename, io="auto")

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

#.........这里部分代码省略.........
    dfi = 2.0 * (nm.pi - open_angle) / nfi
    if is_open:
        nnfi = nfi + 1
    else:
        nnfi = nfi

    is_hollow = force_hollow or not (max(abs(a1), abs(b1)) < 1e-15)

    if is_hollow:
        mr = 0
    else:
        mr = (nnfi - 1) * nl

    grid = nm.zeros((nr, nnfi, nl), dtype=nm.int32)

    n_nod = nr * nnfi * nl - mr
    coors = nm.zeros((n_nod, 3), dtype=nm.float64)

    angles = nm.linspace(open_angle, open_angle+(nfi)*dfi, nfi+1)
    xs = nm.linspace(0.0, length, nl)
    if non_uniform:
        ras = nm.zeros((nr,), dtype=nm.float64)
        rbs = nm.zeros_like(ras)
        advol = (a2**2 - a1**2) / (nr - 1)
        bdvol = (b2**2 - b1**2) / (nr - 1)
        ras[0], rbs[0] = a1, b1
        for ii in range(1, nr):
            ras[ii] = nm.sqrt(advol + ras[ii-1]**2)
            rbs[ii] = nm.sqrt(bdvol + rbs[ii-1]**2)
    else:
        ras = nm.linspace(a1, a2, nr)
        rbs = nm.linspace(b1, b2, nr)

       
##     print dfi * 180.0 / nm.pi
##     print angles * 180.0 / nm.pi
##     print xs
##     print ras
##     print rbs

    # This is 3D only...
    bar = MyBar( "       nodes:" )
    bar.init( n_nod )
    ii = 0
    for ix in range(nr): 
        a, b = ras[ix], rbs[ix]
        for iy, fi in enumerate(angles[:nnfi]):
#            print iy, fi * 180.0 / nm.pi
            for iz, x in enumerate(xs):
##                 print ix, iy, iz, ii
                grid[ix,iy,iz] = ii
                coors[ii] = origin + [x, a * nm.cos(fi), b * nm.sin(fi)]
                if not (ii % 100):
                    bar.update( ii )
                ii += 1

                if not is_hollow and (ix == 0):
                    if iy > 0:
                        grid[ix,iy,iz] = grid[ix,0,iz]
                        ii -= 1
    print
    assert_(ii == n_nod)

    n_el = (nr - 1) * nnfi * (nl - 1)
    conn = nm.zeros((n_el, 8), dtype=nm.int32)

    bar = MyBar( "       elements:" )
    bar.init(n_el)
    ii = 0
    for (ix, iy, iz) in cycle([nr-1, nnfi, nl-1]):
#        print ii, ix, iy, iz
        if iy < (nnfi - 1):
            conn[ii,:] = [grid[ix  ,iy  ,iz  ], grid[ix+1,iy  ,iz  ],
                          grid[ix+1,iy+1,iz  ], grid[ix  ,iy+1,iz  ],
                          grid[ix  ,iy  ,iz+1], grid[ix+1,iy  ,iz+1],
                          grid[ix+1,iy+1,iz+1], grid[ix  ,iy+1,iz+1]]
            ii += 1
        elif not is_open:
            conn[ii,:] = [grid[ix  ,iy  ,iz  ], grid[ix+1,iy  ,iz  ],
                          grid[ix+1,0,iz  ], grid[ix  ,0,iz  ],
                          grid[ix  ,iy  ,iz+1], grid[ix+1,iy  ,iz+1],
                          grid[ix+1,0,iz+1], grid[ix  ,0,iz+1]]
            ii += 1

        if not (ii % 100):
            bar.update( ii )
    print
    mat_id = nm.zeros( (n_el,), dtype = nm.int32 )
    desc = '3_8'

##     print n_nod, n_el, conn.max()
    assert_(n_nod == (conn.max() + 1))

    if axis == 'z':
        coors = coors[:,[1,2,0]]
    elif axis == 'y':
        coors = coors[:,[2,0,1]]

    mesh = Mesh.from_data(name, coors, None, [conn], [mat_id], [desc])
    return mesh

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def gen_block_mesh(dims, shape, centre, name='block'):
    """Generate a 2D or 3D block mesh. The dimension is determined by the
    lenght of the shape argument. 

    Parameters
    ----------

    dims : array of 2 or 3 floats
        Dimensions of the block.
    shape : array of 2 or 3 ints
        Shape (counts of nodes in x, y, z) of the block mesh.
    centre : array of 2 or 3 floats
        Centre of the block.

    name : string
        Mesh name.

    Returns
    -------

    mesh : Mesh instance
    """
    dim = shape.shape[0]

    centre = centre[:dim]
    dims = dims[:dim]

    x0 = centre - 0.5 * dims
    dd = dims / (shape - 1)

    grid = nm.zeros( shape, dtype = nm.int32 )
    n_nod = nm.prod( shape )
    coors = nm.zeros( (n_nod, dim), dtype = nm.float64 )
    
    bar = MyBar( "       nodes:" )
    bar.init( n_nod )
    for ii, ic in enumerate( cycle( shape ) ):
        grid[tuple(ic)] = ii
        coors[ii] = x0 + ic * dd
        if not (ii % 100):
            bar.update( ii )
    bar.update(ii + 1)

    n_el = nm.prod( shape - 1 )
    mat_id = nm.zeros( (n_el,), dtype = nm.int32 )

    if (dim == 2):
        conn = nm.zeros( (n_el, 4), dtype = nm.int32 )
        bar = MyBar( "       elements:" )
        bar.init( n_el )
        for ii, (ix, iy) in enumerate( cycle( shape - 1 ) ):
            conn[ii,:] = [grid[ix  ,iy], grid[ix+1,iy  ],
                          grid[ix+1,iy+1], grid[ix  ,iy+1]]
            if not (ii % 100):
                bar.update( ii )
        bar.update(ii + 1)
        desc = '2_4'

    else:
        conn = nm.zeros( (n_el, 8), dtype = nm.int32 )
        bar = MyBar( "       elements:" )
        bar.init( n_el )
        for ii, (ix, iy, iz) in enumerate( cycle( shape - 1 ) ):
            conn[ii,:] = [grid[ix  ,iy  ,iz  ], grid[ix+1,iy  ,iz  ],
                          grid[ix+1,iy+1,iz  ], grid[ix  ,iy+1,iz  ],
                          grid[ix  ,iy  ,iz+1], grid[ix+1,iy  ,iz+1],
                          grid[ix+1,iy+1,iz+1], grid[ix  ,iy+1,iz+1]]
            if not (ii % 100):
                bar.update( ii )
        bar.update(ii + 1)
        desc = '3_8'

    mesh = Mesh.from_data(name, coors, None, [conn], [mat_id], [desc])
    return mesh

# 需要导入模块: from sfepy.fem import Mesh [as 别名]
# 或者: from sfepy.fem.Mesh import from_data [as 别名]
def refine_reference(geometry, level):
    """
    Refine reference element given by `geometry`.

    Notes
    -----
    The error edges must be generated in the order of the connectivity
    of the previous (lower) level.
    """
    from sfepy.fem import Domain
    from sfepy.fem.geometry_element import geometry_data

    gcoors, gconn = geometry.coors, geometry.conn
    if level == 0:
        return gcoors, gconn

    gd = geometry_data[geometry.name]
    conn = nm.array([gd.conn], dtype=nm.int32)
    mat_id = conn[:, 0].copy()
    mat_id[:] = 0

    mesh = Mesh.from_data("aux", gd.coors, None, [conn], [mat_id], [geometry.name])
    domain = Domain("aux", mesh)

    for ii in range(level):
        domain = domain.refine()

    coors = domain.mesh.coors
    conn = domain.mesh.conns[0]

    n_el = conn.shape[0]

    if geometry.name == "2_3":
        aux_conn = conn.reshape((n_el / 4, 4, 3))

        ir = [[0, 1, 2], [2, 2, 3], [3, 3, 0]]
        ic = [[0, 0, 0], [0, 1, 0], [0, 1, 0]]

    elif geometry.name == "2_4":
        aux_conn = conn.reshape((n_el / 4, 4, 4))

        ir = [[0, 0, 1], [1, 1, 2], [2, 2, 3], [3, 3, 0], [0, 0, 2], [3, 3, 1]]
        ic = [[0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0], [1, 2, 1], [1, 2, 1]]

    elif geometry.name == "3_4":
        aux_conn = conn.reshape((n_el / 8, 8, 4))

        ir = [[0, 0, 1], [1, 1, 2], [2, 0, 0], [3, 1, 1], [3, 2, 2], [3, 0, 0]]
        ic = [[0, 1, 1], [1, 2, 2], [2, 2, 0], [3, 3, 1], [3, 3, 2], [3, 3, 0]]

    elif geometry.name == "3_8":
        aux_conn = conn.reshape((n_el / 8, 8, 8))

        ir = [
            [0, 0, 1],
            [1, 1, 2],
            [2, 2, 3],
            [3, 0, 0],
            [0, 0, 2],
            [0, 0, 1],
            [0, 0, 1],
            [1, 1, 2],
            [2, 2, 3],
            [3, 0, 0],
            [0, 0, 2],
            [0, 0, 1],
            [4, 4, 5],
            [5, 5, 6],
            [6, 6, 7],
            [7, 4, 4],
            [4, 4, 6],
            [4, 4, 5],
            [0, 0, 4],
            [1, 1, 5],
            [2, 2, 6],
            [3, 3, 7],
            [0, 0, 4],
            [1, 1, 5],
            [2, 2, 6],
            [0, 0, 4],
            [0, 0, 4],
        ]
        ic = [
            [0, 1, 0],
            [0, 1, 0],
            [0, 1, 0],
            [0, 3, 0],
            [1, 2, 1],
            [3, 2, 1],
            [4, 5, 4],
            [4, 5, 4],
            [4, 5, 4],
            [4, 7, 4],
            [5, 6, 5],
            [7, 6, 5],
            [0, 3, 0],
            [0, 3, 0],
            [0, 3, 0],
            [0, 1, 0],
            [3, 2, 3],
#.........这里部分代码省略.........

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

#.........这里部分代码省略.........
    o2 = o1 + e_centres.shape[0]
    o3 = o2 + f_centres.shape[0]

    ecc = cmesh.get_conn(cmesh.dim, cmesh.dim - 2)
    eoffs = ecc.offsets

    fcc = cmesh.get_conn(cmesh.dim, cmesh.dim - 1)
    foffs = fcc.offsets

    st = nm.vstack

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        off0, off1 = mesh_in.el_offsets[ig : ig + 2]
        n_el = conn.shape[0]

        e_nodes = ecc.indices[eoffs[off0] : eoffs[off1]].reshape((n_el, 12)) + o1
        f_nodes = fcc.indices[foffs[off0] : foffs[off1]].reshape((n_el, 6)) + o2
        nodes = nm.arange(n_el) + off0 + o3

        c = nm.c_[conn, e_nodes, f_nodes, nodes].T

        new_conn = st(
            [
                c[0],
                c[8],
                c[20],
                c[11],
                c[16],
                c[22],
                c[26],
                c[21],
                c[1],
                c[9],
                c[20],
                c[8],
                c[17],
                c[24],
                c[26],
                c[22],
                c[2],
                c[10],
                c[20],
                c[9],
                c[18],
                c[25],
                c[26],
                c[24],
                c[3],
                c[11],
                c[20],
                c[10],
                c[19],
                c[21],
                c[26],
                c[25],
                c[4],
                c[15],
                c[23],
                c[12],
                c[16],
                c[21],
                c[26],
                c[22],
                c[5],
                c[12],
                c[23],
                c[13],
                c[17],
                c[22],
                c[26],
                c[24],
                c[6],
                c[13],
                c[23],
                c[14],
                c[18],
                c[24],
                c[26],
                c[25],
                c[7],
                c[14],
                c[23],
                c[15],
                c[19],
                c[25],
                c[26],
                c[21],
            ]
        ).T
        new_conn = new_conn.reshape((8 * n_el, 8))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(8)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + "_r", coors, None, conns, mat_ids, mesh_in.descs)

    return mesh