Python mesh_generators.gen_block_mesh函数代码示例

def gen_two_bodies(dims0, shape0, centre0, dims1, shape1, centre1, shift1):
    from sfepy.discrete.fem import Mesh
    from sfepy.mesh.mesh_generators import gen_block_mesh

    m0 = gen_block_mesh(dims0, shape0, centre0)
    m1 = gen_block_mesh(dims1, shape1, centre1)

    coors = nm.concatenate((m0.coors, m1.coors + shift1), axis=0)

    desc = m0.descs[0]
    c0 = m0.get_conn(desc)
    c1 = m1.get_conn(desc)
    conn = nm.concatenate((c0, c1 + m0.n_nod), axis=0)

    ngroups = nm.zeros(coors.shape[0], dtype=nm.int32)
    ngroups[m0.n_nod:] = 1

    mat_id = nm.zeros(conn.shape[0], dtype=nm.int32)
    mat_id[m0.n_el:] = 1

    name = 'two_bodies.mesh'

    mesh = Mesh.from_data(name, coors, ngroups, [conn], [mat_id], m0.descs)

    mesh.write(name, io='auto')

    return mesh

def mesh_hook(mesh, mode):
    if mode == 'read':
        mesh = gen_block_mesh([0.0098, 0.0011, 0.1], [5, 3, 17],
                              [0, 0, 0.05], name='specimen',
                              verbose=False)
        return mesh

    elif mode == 'write':
        pass

    def test_gen_block_mesh(self):
        from sfepy.mesh.mesh_generators import gen_block_mesh

        mesh = gen_block_mesh([1, 2, 3], [4, 3, 5], [2, 1, 0], verbose=False)
        filename = op.join(self.options.out_dir, 'gen_block.mesh')
        mesh.write(filename)

        self.report('block mesh generated')
        return True

    def test_projection_iga_fem(self):
        from sfepy.discrete import FieldVariable
        from sfepy.discrete.fem import FEDomain, Field
        from sfepy.discrete.iga.domain import IGDomain
        from sfepy.mesh.mesh_generators import gen_block_mesh
        from sfepy.discrete.iga.domain_generators import gen_patch_block_domain
        from sfepy.discrete.projections import (make_l2_projection,
                                                make_l2_projection_data)

        shape = [10, 12, 12]
        dims = [5, 6, 6]
        centre = [0, 0, 0]
        degrees = [2, 2, 2]

        nurbs, bmesh, regions = gen_patch_block_domain(dims, shape, centre,
                                                       degrees,
                                                       cp_mode='greville',
                                                       name='iga')
        ig_domain = IGDomain('iga', nurbs, bmesh, regions=regions)

        ig_omega = ig_domain.create_region('Omega', 'all')
        ig_field = Field.from_args('iga', nm.float64, 1, ig_omega,
                                   approx_order='iga', poly_space_base='iga')
        ig_u = FieldVariable('ig_u', 'parameter', ig_field,
                             primary_var_name='(set-to-None)')

        mesh = gen_block_mesh(dims, shape, centre, name='fem')
        fe_domain = FEDomain('fem', mesh)

        fe_omega = fe_domain.create_region('Omega', 'all')
        fe_field = Field.from_args('fem', nm.float64, 1, fe_omega,
                                   approx_order=2)
        fe_u = FieldVariable('fe_u', 'parameter', fe_field,
                             primary_var_name='(set-to-None)')

        def _eval_data(ts, coors, mode, **kwargs):
            return nm.prod(coors**2, axis=1)[:, None, None]

        make_l2_projection_data(ig_u, _eval_data)

        make_l2_projection(fe_u, ig_u) # This calls ig_u.evaluate_at().

        coors = 0.5 * nm.random.rand(20, 3) * dims

        ig_vals = ig_u.evaluate_at(coors)
        fe_vals = fe_u.evaluate_at(coors)

        ok = nm.allclose(ig_vals, fe_vals, rtol=0.0, atol=1e-12)
        if not ok:
            self.report('iga-fem projection failed!')
            self.report('coors:')
            self.report(coors)
            self.report('iga fem diff:')
            self.report(nm.c_[ig_vals, fe_vals, nm.abs(ig_vals - fe_vals)])

        return ok

    def test_gen_block_mesh(self):
        from sfepy.mesh.mesh_generators import gen_block_mesh

        mesh = gen_block_mesh([1, 2, 3], [4, 3, 5], [2, 1, 0], verbose=False)
        filename = op.join(self.options.out_dir, 'gen_block.mesh')
        mesh.write(filename)
        self.report('block mesh generated')

        csum = nm.sum(mesh.coors - nm.min(mesh.coors, axis=0), axis=0)
        return nm.linalg.norm(csum - nm.array([30, 60, 90])) < tolerance

    def mesh_hook(mesh, mode):
        """
        Generate the block mesh.
        """
        if mode == 'read':
            mesh = gen_block_mesh([2, 2, 3], [2, 2, 4], [0, 0, 1.5], name='el3',
                                  verbose=False)
            return mesh

        elif mode == 'write':
            pass

def mesh_hook(mesh, mode):
    """
    Generate the block mesh.
    """
    if mode == 'read':
        mesh = gen_block_mesh(dims, shape, [0, 0], name='user_block',
                              verbose=False)
        return mesh

    elif mode == 'write':
        pass

    def _get_mesh(self, shape):
        """
        Generate an Sfepy rectangular mesh

        Args:
          shape: proposed shape of domain (vertex shape) (n_x, n_y)

        Returns:
          Sfepy mesh

        """
        center = np.zeros_like(shape)
        return gen_block_mesh(shape, np.array(shape) + 1, center,
                              verbose=False)

def main():
    parser = ArgumentParser(description=__doc__,
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('-d', '--dims', metavar='dims',
                        action='store', dest='dims',
                        default='[1.0, 1.0]', help=helps['dims'])
    parser.add_argument('-c', '--centre', metavar='centre',
                        action='store', dest='centre',
                        default='[0.0, 0.0]', help=helps['centre'])
    parser.add_argument('-s', '--shape', metavar='shape',
                        action='store', dest='shape',
                        default='[11, 11]', help=helps['shape'])
    parser.add_argument('--show',
                        action="store_true", dest='show',
                        default=False, help=helps['show'])
    options = parser.parse_args()

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

    output('dimensions:', dims)
    output('centre:    ', centre)
    output('shape:     ', shape)

    mesh = gen_block_mesh(dims, shape, centre, name='block-fem')
    fe_domain = FEDomain('domain', mesh)

    pb, state = run(fe_domain, 1)
    pb.save_state('laplace_shifted_periodic.vtk', state)

    if options.show:
        from sfepy.postprocess.viewer import Viewer
        from sfepy.postprocess.domain_specific import DomainSpecificPlot

        view = Viewer('laplace_shifted_periodic.vtk')
        view(rel_scaling=1,
             domain_specific={'u' : DomainSpecificPlot('plot_warp_scalar',
                                                       ['rel_scaling=1'])},
             is_scalar_bar=True, is_wireframe=True,
             opacity=0.3)

def main():
    parser = OptionParser(usage=usage, version="%prog")
    parser.add_option(
        "-d", "--dims", metavar="dims", action="store", dest="dims", default="[1.0, 1.0]", help=helps["dims"]
    )
    parser.add_option(
        "-c", "--centre", metavar="centre", action="store", dest="centre", default="[0.0, 0.0]", help=helps["centre"]
    )
    parser.add_option(
        "-s", "--shape", metavar="shape", action="store", dest="shape", default="[11, 11]", help=helps["shape"]
    )
    parser.add_option("", "--show", action="store_true", dest="show", default=False, help=helps["show"])
    (options, args) = parser.parse_args()

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

    output("dimensions:", dims)
    output("centre:    ", centre)
    output("shape:     ", shape)

    mesh = gen_block_mesh(dims, shape, centre, name="block-fem")
    fe_domain = FEDomain("domain", mesh)

    pb, state = run(fe_domain, 1)
    pb.save_state("laplace_shifted_periodic.vtk", state)

    if options.show:
        from sfepy.postprocess.viewer import Viewer
        from sfepy.postprocess.domain_specific import DomainSpecificPlot

        view = Viewer("laplace_shifted_periodic.vtk")
        view(
            rel_scaling=1,
            domain_specific={"u": DomainSpecificPlot("plot_warp_scalar", ["rel_scaling=1"])},
            is_scalar_bar=True,
            is_wireframe=True,
            opacity=0.3,
        )

def main():
    parser = ArgumentParser(description=__doc__)
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('-o', metavar='filename',
                        action='store', dest='output_filename',
                        default='out.vtk', help=helps['filename'])
    parser.add_argument('-f', '--format', metavar='format',
                        action='store', type=str, dest='format',
                        default=None, help=helps['format'])
    parser.add_argument('-d', '--dims', metavar='dims',
                        action='store', dest='dims',
                        default='[1.0, 1.0, 1.0]', help=helps['dims'])
    parser.add_argument('-s', '--shape', metavar='shape',
                        action='store', dest='shape',
                        default='[11, 11, 11]', help=helps['shape'])
    parser.add_argument('-c', '--centre', metavar='centre',
                        action='store', dest='centre',
                        default='[0.0, 0.0, 0.0]', help=helps['centre'])
    parser.add_argument('-2', '--2d',
                        action='store_true', dest='is_2d',
                        default=False, help=helps['2d'])
    options = parser.parse_args()

    dim = 2 if options.is_2d else 3

    dims = nm.array(eval(options.dims), dtype=nm.float64)[:dim]
    shape = nm.array(eval(options.shape), dtype=nm.int32)[:dim]
    centre = nm.array(eval(options.centre), dtype=nm.float64)[:dim]

    output.prefix = 'blockgen:'
    output('dimensions:', dims)
    output('shape:', shape)
    output('centre:', centre)

    mesh = gen_block_mesh(dims, shape, centre, name=options.output_filename)

    io = MeshIO.for_format(options.output_filename, format=options.format,
                           writable=True)

    mesh.write(options.output_filename, io=io)

def make_mesh(dims, shape, transform=None):
    """
    Generate a 2D rectangle mesh in 3D space, and optionally apply a coordinate
    transform.
    """
    _mesh = gen_block_mesh(dims, shape, [0, 0], name='shell10x', verbose=False)

    coors = nm.c_[_mesh.coors, nm.zeros(_mesh.n_nod, dtype=nm.float64)]
    coors = nm.ascontiguousarray(coors)

    conns = [_mesh.get_conn(_mesh.descs[0])]

    mesh = Mesh.from_data(_mesh.name, coors, _mesh.cmesh.vertex_groups, conns,
                          [_mesh.cmesh.cell_groups], _mesh.descs)

    if transform == 'bend':
        bbox = mesh.get_bounding_box()
        x0, x1 = bbox[:, 0]

        angles = 0.5 *  nm.pi * (coors[:, 0] - x0) / (x1 - x0)
        mtx = make_axis_rotation_matrix([0, -1, 0], angles[:, None, None])

        coors = mesh.coors.copy()
        coors[:, 0] = 0
        coors[:, 2] = (x1 - x0)

        mesh.coors[:] = transform_data(coors, mtx=mtx)
        mesh.coors[:, 0] -= 0.5 * (x1 - x0)

    elif transform == 'twist':
        bbox = mesh.get_bounding_box()
        x0, x1 = bbox[:, 0]

        angles = 0.5 *  nm.pi * (coors[:, 0] - x0) / (x1 - x0)
        mtx = make_axis_rotation_matrix([-1, 0, 0], angles[:, None, None])

        mesh.coors[:] = transform_data(mesh.coors, mtx=mtx)

    return mesh

    def test_laplace_shifted_periodic(self):
        import numpy as nm
        from sfepy.mesh.mesh_generators import gen_block_mesh
        from sfepy.discrete.fem import FEDomain
        from examples.diffusion.laplace_shifted_periodic import run

        dims = [2.0, 1.0]
        shape = [21, 11]
        centre = [0.0, 0.0]
        mesh = gen_block_mesh(dims, shape, centre, name='block-fem')
        fe_domain = FEDomain('domain', mesh)

        pb, state = run(fe_domain, 3)

        gamma3 = pb.domain.regions['Gamma3']
        gamma4 = pb.domain.regions['Gamma4']

        field = pb.fields['fu']

        # Check that the shift equals to one.
        i3 = field.get_dofs_in_region(gamma3, merge=True)
        i4 = field.get_dofs_in_region(gamma4, merge=True)

        i_corners = nm.array([0, shape[0] - 1])
        ii = nm.setdiff1d(nm.arange(len(i3)), i_corners)

        vals = state()

        shift = vals[i3] - vals[i4]

        ok = (shift[i_corners] == 0.0).all()

        ok = ok and nm.allclose(shift[ii], 1.0, rtol=0.0, atol=1e-14)

        if not ok:
            self.report('wrong shift:', shift)

        return ok

def main():
    parser = ArgumentParser(description=__doc__.rstrip(),
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('output_dir', help=helps['output_dir'])
    parser.add_argument('--dims', metavar='dims',
                        action='store', dest='dims',
                        default='1.0,1.0,1.0', help=helps['dims'])
    parser.add_argument('--shape', metavar='shape',
                        action='store', dest='shape',
                        default='7,7,7', help=helps['shape'])
    parser.add_argument('--centre', metavar='centre',
                        action='store', dest='centre',
                        default='0.0,0.0,0.0', help=helps['centre'])
    parser.add_argument('-3', '--3d',
                        action='store_true', dest='is_3d',
                        default=False, help=helps['3d'])
    parser.add_argument('--order', metavar='int', type=int,
                        action='store', dest='order',
                        default=1, help=helps['order'])
    options = parser.parse_args()

    dim = 3 if options.is_3d else 2
    dims = nm.array(eval(options.dims), dtype=nm.float64)[:dim]
    shape = nm.array(eval(options.shape), dtype=nm.int32)[:dim]
    centre = nm.array(eval(options.centre), dtype=nm.float64)[:dim]

    output('dimensions:', dims)
    output('shape:     ', shape)
    output('centre:    ', centre)

    mesh0 = gen_block_mesh(dims, shape, centre, name='block-fem',
                           verbose=True)
    domain0 = FEDomain('d', mesh0)

    bbox = domain0.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    eps = 1e-8 * (max_x - min_x)

    cnt = (shape[0] - 1) // 2
    g0 = 0.5 * dims[0]
    grading = nm.array([g0 / 2**ii for ii in range(cnt)]) + eps + centre[0] - g0

    domain, subs = refine_towards_facet(domain0, grading, 'x <')

    omega = domain.create_region('Omega', 'all')

    gamma1 = domain.create_region('Gamma1',
                                  'vertices in (x < %.10f)' % (min_x + eps),
                                  'facet')
    gamma2 = domain.create_region('Gamma2',
                                  'vertices in (x > %.10f)' % (max_x - eps),
                                  'facet')

    field = Field.from_args('fu', nm.float64, 1, omega,
                            approx_order=options.order)

    if subs is not None:
        field.substitute_dofs(subs)

    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')

    integral = Integral('i', order=2*options.order)

    t1 = Term.new('dw_laplace(v, u)',
                  integral, omega, v=v, u=u)
    eq = Equation('eq', t1)
    eqs = Equations([eq])

    def u_fun(ts, coors, bc=None, problem=None):
        """
        Define a displacement depending on the y coordinate.
        """
        if coors.shape[1] == 2:
            min_y, max_y = bbox[:, 1]
            y = (coors[:, 1] - min_y) / (max_y - min_y)

            val = (max_y - min_y) * nm.cos(3 * nm.pi * y)

        else:
            min_y, max_y = bbox[:, 1]
            min_z, max_z = bbox[:, 2]
            y = (coors[:, 1] - min_y) / (max_y - min_y)
            z = (coors[:, 2] - min_z) / (max_z - min_z)

            val = ((max_y - min_y) * (max_z - min_z)
                   * nm.cos(3 * nm.pi * y) * (1.0 + 3.0 * (z - 0.5)**2))

        return val

    bc_fun = Function('u_fun', u_fun)
    fix1 = EssentialBC('shift_u', gamma1, {'u.0' : bc_fun})
    fix2 = EssentialBC('fix2', gamma2, {'u.all' : 0.0})

    ls = ScipyDirect({})

    nls = Newton({}, lin_solver=ls)

    pb = Problem('heat', equations=eqs, nls=nls, ls=ls)

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

def main():
    parser = ArgumentParser(description=__doc__.rstrip(),
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('output_dir', help=helps['output_dir'])
    parser.add_argument('--dims', metavar='dims',
                        action='store', dest='dims',
                        default='1.0,1.0,1.0', help=helps['dims'])
    parser.add_argument('--shape', metavar='shape',
                        action='store', dest='shape',
                        default='11,11,11', help=helps['shape'])
    parser.add_argument('--centre', metavar='centre',
                        action='store', dest='centre',
                        default='0.0,0.0,0.0', help=helps['centre'])
    parser.add_argument('-2', '--2d',
                        action='store_true', dest='is_2d',
                        default=False, help=helps['2d'])
    parser.add_argument('--u-order', metavar='int', type=int,
                        action='store', dest='order_u',
                        default=1, help=helps['u-order'])
    parser.add_argument('--p-order', metavar='int', type=int,
                        action='store', dest='order_p',
                        default=1, help=helps['p-order'])
    parser.add_argument('--linearization', choices=['strip', 'adaptive'],
                        action='store', dest='linearization',
                        default='strip', help=helps['linearization'])
    parser.add_argument('--metis',
                        action='store_true', dest='metis',
                        default=False, help=helps['metis'])
    parser.add_argument('--silent',
                        action='store_true', dest='silent',
                        default=False, help=helps['silent'])
    parser.add_argument('--clear',
                        action='store_true', dest='clear',
                        default=False, help=helps['clear'])
    options, petsc_opts = parser.parse_known_args()

    comm = pl.PETSc.COMM_WORLD

    output_dir = options.output_dir

    filename = os.path.join(output_dir, 'output_log_%02d.txt' % comm.rank)
    if comm.rank == 0:
        ensure_path(filename)
    comm.barrier()

    output.prefix = 'sfepy_%02d:' % comm.rank
    output.set_output(filename=filename, combined=options.silent == False)

    output('petsc options:', petsc_opts)

    mesh_filename = os.path.join(options.output_dir, 'para.h5')

    if comm.rank == 0:
        from sfepy.mesh.mesh_generators import gen_block_mesh

        if options.clear:
            for _f in chain(*[glob.glob(os.path.join(output_dir, clean_pattern))
                              for clean_pattern
                              in ['*.h5', '*.txt', '*.png']]):
                output('removing "%s"' % _f)
                os.remove(_f)

        dim = 2 if options.is_2d else 3
        dims = nm.array(eval(options.dims), dtype=nm.float64)[:dim]
        shape = nm.array(eval(options.shape), dtype=nm.int32)[:dim]
        centre = nm.array(eval(options.centre), dtype=nm.float64)[:dim]

        output('dimensions:', dims)
        output('shape:     ', shape)
        output('centre:    ', centre)

        mesh = gen_block_mesh(dims, shape, centre, name='block-fem',
                              verbose=True)
        mesh.write(mesh_filename, io='auto')

    comm.barrier()

    output('field u order:', options.order_u)
    output('field p order:', options.order_p)

    solve_problem(mesh_filename, options, comm)