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)

def main(cli_args):
    dims = parse_argument_list(cli_args.dims, float)
    shape = parse_argument_list(cli_args.shape, int)
    centre = parse_argument_list(cli_args.centre, float)
    material_parameters = parse_argument_list(cli_args.material_parameters,
                                              float)
    order = cli_args.order

    ts_vals = cli_args.ts.split(',')
    ts = {
        't0' : float(ts_vals[0]), 't1' : float(ts_vals[1]),
        'n_step' : int(ts_vals[2])}

    do_plot = cli_args.plot

    ### Mesh and regions ###
    mesh = gen_block_mesh(
        dims, shape, centre, name='block', verbose=False)
    domain = FEDomain('domain', mesh)

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

    lbn, rtf = domain.get_mesh_bounding_box()
    box_regions = define_box_regions(3, lbn, rtf)
    regions = dict([
        [r, domain.create_region(r, box_regions[r][0], box_regions[r][1])]
        for r in box_regions])

    ### Fields ###
    scalar_field = Field.from_args(
        'fu', np.float64, 'scalar', omega, approx_order=order-1)
    vector_field = Field.from_args(
        'fv', np.float64, 'vector', omega, approx_order=order)

    u = FieldVariable('u', 'unknown', vector_field, history=1)
    v = FieldVariable('v', 'test', vector_field, primary_var_name='u')
    p = FieldVariable('p', 'unknown', scalar_field, history=1)
    q = FieldVariable('q', 'test', scalar_field, primary_var_name='p')

    ### Material ###
    c10, c01 = material_parameters
    m = Material(
        'm', mu=2*c10, kappa=2*c01,
    )

    ### Boundary conditions ###
    x_sym = EssentialBC('x_sym', regions['Left'], {'u.0' : 0.0})
    y_sym = EssentialBC('y_sym', regions['Near'], {'u.1' : 0.0})
    z_sym = EssentialBC('z_sym', regions['Bottom'], {'u.2' : 0.0})
    disp_fun = Function('disp_fun', get_displacement)
    displacement = EssentialBC(
        'displacement', regions['Right'], {'u.0' : disp_fun})
    ebcs = Conditions([x_sym, y_sym, z_sym, displacement])

    ### Terms and equations ###
    integral = Integral('i', order=2*order)

    term_neohook = Term.new(
        'dw_tl_he_neohook(m.mu, v, u)',
        integral, omega, m=m, v=v, u=u)
    term_mooney = Term.new(
        'dw_tl_he_mooney_rivlin(m.kappa, v, u)',
        integral, omega, m=m, v=v, u=u)
    term_pressure = Term.new(
        'dw_tl_bulk_pressure(v, u, p)',
        integral, omega, v=v, u=u, p=p)

    term_volume_change = Term.new(
        'dw_tl_volume(q, u)',
        integral, omega, q=q, u=u, term_mode='volume')
    term_volume = Term.new(
        'dw_volume_integrate(q)',
        integral, omega, q=q)

    eq_balance = Equation('balance', term_neohook+term_mooney+term_pressure)
    eq_volume = Equation('volume', term_volume_change-term_volume)
    equations = Equations([eq_balance, eq_volume])

    ### Solvers ###
    ls = ScipyDirect({})
    nls_status = IndexedStruct()
    nls = Newton(
        {'i_max' : 5},
        lin_solver=ls, status=nls_status
    )

    ### Problem ###
    pb = Problem('hyper', equations=equations)
    pb.set_bcs(ebcs=ebcs)
    pb.set_ics(ics=Conditions([]))
    tss = SimpleTimeSteppingSolver(ts, nls=nls, context=pb)
    pb.set_solver(tss)

    ### Solution ###
    axial_stress = []
    axial_displacement = []
    def stress_strain_fun(*args, **kwargs):
        return stress_strain(
            *args, order=order, global_stress=axial_stress,
            global_displacement=axial_displacement, **kwargs)
#.........这里部分代码省略.........

    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

import scipy

aux = "eig.scipy,method:'eigh',tol:1e-7,maxiter:10000".split(',')
kwargs = {}
for option in aux[1:]:
    key, val = option.split(':')
    kwargs[key.strip()] = eval(val)
eig_conf = Struct(name='evp', kind=aux[0], **kwargs)

dims = nm.array([1.0, 1.5], dtype=nm.float64)
dim = len(dims)

centre = nm.array([0.0, 0.0], dtype=nm.float64)[:dim]
shape = nm.array([11, 16], dtype=nm.int32)[:dim]

mesh = gen_block_mesh(dims, shape, centre, name='mesh')

eig_solver = Solver.any_from_conf(eig_conf)

# Build the problem definition.
domain = FEDomain('domain', mesh)

bbox = domain.get_mesh_bounding_box()
min_coor, max_coor = bbox[:, -1]
eps = 1e-8 * (max_coor - min_coor)
ax = 'xyz'[:dim][-1]

omega = domain.create_region('Omega', 'all')
bottom = domain.create_region('Bottom', 'vertices in (%s < %.10f)' % (ax, min_coor + eps), 'facet')
bottom_top = domain.create_region('BottomTop', 'r.Bottom +v vertices in (%s > %.10f)' % (ax, max_coor - eps), 'facet')

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)

def main():
    from sfepy import data_dir

    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-s', '--show',
                      action="store_true", dest='show',
                      default=False, help=help['show'])
    options, args = parser.parse_args()
    options_probe = True
    folder = str(uuid.uuid4())
    os.mkdir(folder)
    os.chdir(folder)
    
    file = open('README.txt', 'w')
    file.write('DIMENSIONS\n')
    file.write('Lx = '+str(dims[0])+' Ly = '+str(dims[1])+' Lz = '+str(dims[2])+'\n')
    file.write('DISCRETIZATION (NX, NY, NZ)\n')
    file.write(str(NX)+'  '+str(NY)+'  '+str(NZ)+'\n')
    file.write('MATERIALS\n')
    file.write(str(E_f)+' '+str(nu_f)+' '+str(E_m)+' '+str(nu_m)+'\n')
    
    #mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh')
    
    mesh = mesh_generators.gen_block_mesh(dims,shape,centre,name='block')
    domain = FEDomain('domain', mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:,0]
    min_y, max_y = domain.get_mesh_bounding_box()[:,1]
    min_z, max_z = domain.get_mesh_bounding_box()[:,2]
    eps = 1e-8 * (max_x - min_x)
    print min_x, max_x
    print min_y, max_y
    print min_z, max_z
    R1 = domain.create_region('Ym',
                                  'vertices in z < %.10f' % (max_z/2))
    R2 = domain.create_region('Yf',
                                  'vertices in z >= %.10f' % (min_z/2))
    omega = domain.create_region('Omega', 'all')
    gamma1 = domain.create_region('Left',
                                  'vertices in x < %.10f' % (min_x + eps), 
                                  'facet')
    gamma2 = domain.create_region('Right',
                                  'vertices in x > %.10f' % (max_x - eps),
                                  'facet')
    gamma3 = domain.create_region('Front',
                                  'vertices in y < %.10f' % (min_y + eps),
                                  'facet')
    gamma4 = domain.create_region('Back',
                                  'vertices in y > %.10f' % (max_y - eps),
                                  'facet')
    gamma5 = domain.create_region('Bottom',
                                  'vertices in z < %.10f' % (min_z + eps),
                                  'facet')
    gamma6 = domain.create_region('Top',
                                  'vertices in z > %.10f' % (max_z - eps),
                                  'facet')



    field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2)

    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')
    mu=1.1
    lam=1.0
    m = Material('m', lam=lam, mu=mu)
    f = Material('f', val=[[0.0], [0.0],[0.0]])
    #mu,lam=m.get_constants_mu_lam()
    #print mu.lam 
    D = stiffness_from_lame(3,lam, mu)    
    mat = Material('Mat', D=D)

    #D = stiffness_from_youngpoisson(2, options.young, options.poisson)
    get_mat = Function('get_mat1',get_mat1)
    #get_mat1=Function('get_mat', (lambda ts, coors, mode=None, problem=None, **kwargs:
    #                get_mat(coors, mode, problem)))
    #mat = Material('Mat', function=Function('get_mat1',get_mat1))
    #mat = Material('Mat', 'get_mat')
    integral = Integral('i', order=3)

    t1 = Term.new('dw_lin_elastic(Mat.D, v, u)',
         integral, omega, Mat=mat, v=v, u=u)
    t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
    eq = Equation('balance', t1 + t2)
    eqs = Equations([eq])

    fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})
    left_bc  = EssentialBC('Left',  gamma1, {'u.0' : 0.0})
    right_bc = EssentialBC('Right', gamma2, {'u.0' : 0.0})
    back_bc = EssentialBC('Front', gamma3, {'u.1' : 0.0})
    front_bc = EssentialBC('Back', gamma4, {'u.1' : 0.0})
    bottom_bc = EssentialBC('Bottom', gamma5, {'u.all' : 0.0})
    top_bc = EssentialBC('Top', gamma6, {'u.2' : 0.2})

    bc=[left_bc,right_bc,back_bc,front_bc,bottom_bc,top_bc]
    #bc=[bottom_bc,top_bc]

    bc_fun = Function('shift_u_fun', shift_u_fun, extra_args={'shift' : 0.01})
    shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})
    #get_mat = Function('get_mat1',get_mat1)
#.........这里部分代码省略.........