def evalValAndDeriv(D):
    m = Material('m', D = D, rho = 2700.0)

    integral = Integral('i', order=2)

    t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=m, v=v, u=u)
    eq1 = Equation('stiffness', t1)
    eq2 = Equation('mass', t2)
    lhs_eqs = Equations([eq1, eq2])

    pb = Problem('modal', equations = lhs_eqs)

    pb.time_update()
    n_rbm = dim * (dim + 1) / 2

    pb.update_materials()

    # Assemble stiffness and mass matrices.
    mtx_k = eq1.evaluate(mode='weak', dw_mode='matrix', asm_obj=pb.mtx_a)
    mtx_m = mtx_k.copy()
    mtx_m.data[:] = 0.0
    mtx_m = eq2.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx_m)

    eigs0, evecs0 = scipy.sparse.linalg.eigsh(mtx_k, k = 10, M = mtx_m, which = 'SM')

    eigs = eigs0[3:]
    evecs = evecs0[:, 3:]

    dydmu = numpy.array([evecs[:, i].T.dot(dKdmu.dot(evecs[:, i])) for i in range(evecs.shape[1])])
    dydlambda = numpy.array([evecs[:, i].T.dot(dKdlambda.dot(evecs[:, i])) for i in range(evecs.shape[1])])

    return eigs, dydmu, dydlambda

def assemble(mtx_d):
    m = Material('m', D=mtx_d, rho=density)

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

    t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=m, v=v, u=u)
    eq1 = Equation('stiffness', t1)
    eq2 = Equation('mass', t2)
    lhs_eqs = Equations([eq1, eq2])

    pb = Problem('modal', equations=lhs_eqs)

    pb.time_update()
    n_rbm = dim * (dim + 1) / 2

    pb.update_materials()

    tmp = time.time()
    # Assemble stiffness and mass matrices.
    mtx_k = eq1.evaluate(mode='weak', dw_mode='matrix', asm_obj=pb.mtx_a)
    mtx_m = mtx_k.copy()
    mtx_m.data[:] = 0.0
    mtx_m = eq2.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx_m)

    return mtx_k, mtx_m

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()

    mesh = Mesh.from_file(data_dir + "/meshes/2d/rectangle_tri.mesh")
    domain = Domain("domain", mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
    eps = 1e-8 * (max_x - min_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, "vector", omega, approx_order=2)

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

    m = Material("m", lam=1.0, mu=1.0)
    f = Material("f", val=[[0.02], [0.01]])

    integral = Integral("i", order=3)

    t1 = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v, u)", integral, omega, m=m, 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})

    bc_fun = Function("shift_u_fun", shift_u_fun, extra_args={"shift": 0.01})
    shift_u = EssentialBC("shift_u", gamma2, {"u.0": bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem("elasticity", equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups("regions")

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print nls_status

    pb.save_state("linear_elasticity.vtk", vec)

    if options.show:
        view = Viewer("linear_elasticity.vtk")
        view(vector_mode="warp_norm", rel_scaling=2, is_scalar_bar=True, is_wireframe=True)

def solve_problem(shape, dims, young, poisson, force, transform=None):
    domain = make_domain(dims[:2], shape, transform=transform)

    omega = domain.regions['Omega']
    gamma1 = domain.regions['Gamma1']
    gamma2 = domain.regions['Gamma2']

    field = Field.from_args('fu', nm.float64, 6, omega, approx_order=1,
                            poly_space_base='shell10x')
    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')

    thickness = dims[2]
    if transform is None:
        pload = [[0.0, 0.0, force / shape[1], 0.0, 0.0, 0.0]] * shape[1]

    elif transform == 'bend':
        pload = [[force / shape[1], 0.0, 0.0, 0.0, 0.0, 0.0]] * shape[1]

    elif transform == 'twist':
        pload = [[0.0, force / shape[1], 0.0, 0.0, 0.0, 0.0]] * shape[1]

    m = Material('m', D=sh.create_elastic_tensor(young=young, poisson=poisson),
                 values={'.drill' : 1e-7})
    load = Material('load', values={'.val' : pload})

    aux = Integral('i', order=3)
    qp_coors, qp_weights = aux.get_qp('3_8')
    qp_coors[:, 2] = thickness * (qp_coors[:, 2] - 0.5)
    qp_weights *= thickness

    integral = Integral('i', coors=qp_coors, weights=qp_weights, order='custom')

    t1 = Term.new('dw_shell10x(m.D, m.drill, v, u)',
                  integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_point_load(load.val, v)',
                  integral, gamma2, load=load, v=v)
    eq = Equation('balance', t1 - t2)
    eqs = Equations([eq])

    fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem('elasticity with shell10x', equations=eqs)
    pb.set_bcs(ebcs=Conditions([fix_u]))
    pb.set_solver(nls)

    state = pb.solve()

    return pb, state, u, gamma2

    def test_invariance_qp(self):
        from sfepy import data_dir
        from sfepy.discrete import Integral
        from sfepy.terms import Term
        from sfepy.discrete.common.mappings import get_physical_qps

        ok = True
        for name in [
            "meshes/3d/block.mesh",
            "meshes/3d/cylinder.mesh",
            "meshes/2d/square_quad.mesh",
            "meshes/2d/square_unit_tri.mesh",
        ]:
            self.report(name)

            u = prepare_variable(op.join(data_dir, name), n_components=3)
            omega = u.field.region

            integral = Integral("i", order=3)
            qps = get_physical_qps(omega, integral)
            coors = qps.values

            term = Term.new("ev_volume_integrate(u)", integral, omega, u=u)
            term.setup()
            val1 = term.evaluate(mode="qp")
            val1 = val1.ravel()

            val2 = u.evaluate_at(coors).ravel()

            self.report("value: max. difference:", nm.abs(val1 - val2).max())
            ok1 = nm.allclose(val1, val2, rtol=0.0, atol=1e-12)
            self.report("->", ok1)

            term = Term.new("ev_grad(u)", integral, omega, u=u)
            term.setup()
            val1 = term.evaluate(mode="qp")
            val1 = val1.ravel()

            val2 = u.evaluate_at(coors, mode="grad").ravel()

            self.report("gradient: max. difference:", nm.abs(val1 - val2).max())
            ok2 = nm.allclose(val1, val2, rtol=0.0, atol=1e-10)
            self.report("->", ok2)

            ok = ok and ok1 and ok2

        return ok

def run(domain, order):
    omega = domain.create_region('Omega', 'all')
    bbox = domain.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    min_y, max_y = bbox[:, 1]
    eps = 1e-8 * (max_x - min_x)
    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')
    gamma3 = domain.create_region('Gamma3',
                                  'vertices in y < %.10f' % (min_y + eps),
                                  'facet')
    gamma4 = domain.create_region('Gamma4',
                                  'vertices in y > %.10f' % (max_y - eps),
                                  'facet')

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

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

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

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

    fix1 = EssentialBC('fix1', gamma1, {'u.0' : 0.4})
    fix2 = EssentialBC('fix2', gamma2, {'u.0' : 0.0})

    def get_shift(ts, coors, region):
        return nm.ones_like(coors[:, 0])

    dof_map_fun = Function('dof_map_fun', per.match_x_line)
    shift_fun = Function('shift_fun', get_shift)

    sper = LinearCombinationBC('sper', [gamma3, gamma4], {'u.0' : 'u.0'},
                               dof_map_fun, 'shifted_periodic',
                               arguments=(shift_fun,))

    ls = ScipyDirect({})

    pb = Problem('laplace', equations=eqs, auto_solvers=None)

    pb.time_update(ebcs=Conditions([fix1, fix2]), lcbcs=Conditions([sper]))

    ev = pb.get_evaluator()
    nls = Newton({}, lin_solver=ls,
                 fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

    pb.set_solver(nls)

    state = pb.solve()

    return pb, state

    def test_save_ebc(self):
        from sfepy.discrete import (FieldVariable, Integral,
                                    Equation, Equations, Problem)
        from sfepy.discrete.conditions import Conditions, EssentialBC
        from sfepy.terms import Term

        name = op.join(self.options.out_dir,
                       op.splitext(op.basename(__file__))[0])

        integral = Integral('i', order=1)

        u = self.variables['u']
        v = FieldVariable('v', 'test', u.field, primary_var_name='u')

        p = self.variables['p']
        q = FieldVariable('q', 'test', p.field, primary_var_name='p')

        regions = self.problem.domain.regions
        omega = regions['Omega']

        # Problem.save_ebc() requires to have equations defined.
        t1 = Term.new('dw_lin_elastic(v, u)',
                      integral, omega, v=v, u=u)
        t2 = Term.new('dw_laplace(q, p)', integral, omega, q=q, p=p)
        eq = Equation('aux', t1 + t2)
        eqs = Equations([eq])

        pb = Problem('test', equations=eqs, auto_solvers=False)

        all_ebcs = []
        all_ebcs.append(EssentialBC('fix_u1', regions['RightFix'],
                                    {'u.all' : nm.array([0.0, 1.0])}))
        all_ebcs.append(EssentialBC('fix_u2', regions['LeftStrip'],
                                    {'u.0' : 0.0, 'u.1' : 1.0}))
        all_ebcs.append(EssentialBC('fix_p1', regions['LeftFix'],
                                    {'p.all' : 0.0}))
        all_ebcs.append(EssentialBC('fix_p2', regions['RightStrip'],
                                    {'p.0' : 0.0}))
        ebcs = Conditions(all_ebcs)

        pb.time_update(ebcs=ebcs)

        pb.save_ebc(name + '_ebcs_f.vtk', ebcs=ebcs, force=True)
        pb.save_ebc(name + '_ebcs.vtk', ebcs=ebcs, default=-1, force=False)

        return True

def make_l2_projection(target, source):
    """
    Project `source` field variable to `target` field variable using
    :math:`L^2` dot product.
    """
    order = target.field.get_true_order()**2
    integral = Integral('i', order=order)

    un = target.name
    v = FieldVariable('v', 'test', target.field, 1, primary_var_name=un)
    lhs = Term.new('dw_mass_scalar(v, %s)' % un, integral,
                   target.field.region, v=v, **{un : target})

    def eval_variable(ts, coors, mode, **kwargs):
        if mode == 'qp':
            val = source.evaluate_at(coors)
            val.shape = val.shape + (1,)
            out = {'val' : val}
            return out

    m = Material('m', function=eval_variable)
    rhs = Term.new('dw_volume_lvf(m.val, v)', integral, target.field.region,
                   m=m, v=v)

    eq = Equation('projection', lhs - rhs)
    eqs = Equations([eq])

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = ProblemDefinition('aux', equations=eqs, nls=nls, ls=ls)

    pb.time_update()

    # This sets the target variable with the projection solution.
    pb.solve()

    if nls_status.condition != 0:
        output('L2 projection: solver did not converge!')

def make_l2_projection_data(target, eval_data):
    """
    Project scalar data given by a material-like `eval_data()` function to a
    scalar `target` field variable using the :math:`L^2` dot product.
    """
    order = target.field.approx_order * 2
    integral = Integral("i", order=order)

    un = target.name
    v = FieldVariable("v", "test", target.field, primary_var_name=un)
    lhs = Term.new("dw_volume_dot(v, %s)" % un, integral, target.field.region, v=v, **{un: target})

    def _eval_data(ts, coors, mode, **kwargs):
        if mode == "qp":
            val = eval_data(ts, coors, mode, **kwargs)
            return {"val": val}

    m = Material("m", function=_eval_data)
    rhs = Term.new("dw_volume_lvf(m.val, v)", integral, target.field.region, m=m, v=v)

    eq = Equation("projection", lhs - rhs)
    eqs = Equations([eq])

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = ProblemDefinition("aux", equations=eqs, nls=nls, ls=ls)

    pb.time_update()

    # This sets the target variable with the projection solution.
    pb.solve()

    if nls_status.condition != 0:
        output("L2 projection: solver did not converge!")

def run(domain, order):
    omega = domain.create_region("Omega", "all")
    bbox = domain.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    min_y, max_y = bbox[:, 1]
    eps = 1e-8 * (max_x - min_x)
    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")
    gamma3 = domain.create_region("Gamma3", "vertices in y < %.10f" % (min_y + eps), "facet")
    gamma4 = domain.create_region("Gamma4", "vertices in y > %.10f" % (max_y - eps), "facet")

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

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

    integral = Integral("i", order=2 * order)

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

    fix1 = EssentialBC("fix1", gamma1, {"u.0": 0.4})
    fix2 = EssentialBC("fix2", gamma2, {"u.0": 0.0})

    def get_shift(ts, coors, region):
        return nm.ones_like(coors[:, 0])

    dof_map_fun = Function("dof_map_fun", per.match_x_line)
    shift_fun = Function("shift_fun", get_shift)

    sper = LinearCombinationBC(
        "sper", [gamma3, gamma4], {"u.0": "u.0"}, dof_map_fun, "shifted_periodic", arguments=(shift_fun,)
    )

    ls = ScipyDirect({})

    pb = Problem("laplace", equations=eqs, auto_solvers=None)

    pb.time_update(ebcs=Conditions([fix1, fix2]), lcbcs=Conditions([sper]))

    ev = pb.get_evaluator()
    nls = Newton({}, lin_solver=ls, fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

    pb.set_solver(nls)

    state = pb.solve()

    return pb, state

    def test_invariance_qp(self):
        from sfepy import data_dir
        from sfepy.discrete import Variables, Integral
        from sfepy.discrete.fem import Mesh, FEDomain, Field
        from sfepy.terms import Term
        from sfepy.discrete.common.mappings import get_physical_qps

        mesh = Mesh.from_file(data_dir + '/meshes/3d/block.mesh')

        bbox = mesh.get_bounding_box()
        dd = bbox[1,:] - bbox[0,:]
        data = nm.sin(4.0 * nm.pi * mesh.coors[:,0:1] / dd[0]) \
               * nm.cos(4.0 * nm.pi * mesh.coors[:,1:2] / dd[1])

        variables = {
            'u'       : ('unknown field', 'scalar_tp', 0),
            'v'       : ('test field',    'scalar_tp', 'u'),
        }

        domain = FEDomain('domain', mesh)
        omega = domain.create_region('Omega', 'all')
        field = Field.from_args('scalar_tp', nm.float64, 1, omega,
                                approx_order=1)
        ff = {field.name : field}

        vv = Variables.from_conf(transform_variables(variables), ff)
        u = vv['u']
        u.set_from_mesh_vertices(data)

        integral = Integral('i', order=2)
        term = Term.new('ev_volume_integrate(u)', integral, omega, u=u)
        term.setup()
        val1 = term.evaluate(mode='qp')
        val1 = val1.ravel()

        qps = get_physical_qps(omega, integral)
        coors = qps.values

        val2 = u.evaluate_at(coors).ravel()

        self.report('max. difference:', nm.abs(val1 - val2).max())
        ok = nm.allclose(val1, val2, rtol=0.0, atol=1e-12)
        self.report('invariance in qp: %s' % ok)

        return ok

def main():
    from sfepy import data_dir

    parser = ArgumentParser()
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('-s', '--show',
                        action="store_true", dest='show',
                        default=False, help=helps['show'])
    options = parser.parse_args()

    mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh')
    domain = FEDomain('domain', mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:,0]
    eps = 1e-8 * (max_x - min_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, 'vector', omega,
                            approx_order=2)

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

    m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))
    f = Material('f', val=[[0.02], [0.01]])

    integral = Integral('i', order=3)

    t1 = Term.new('dw_lin_elastic(m.D, v, u)',
                  integral, omega, m=m, 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})

    bc_fun = Function('shift_u_fun', shift_u_fun,
                      extra_args={'shift' : 0.01})
    shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem('elasticity', equations=eqs)
    pb.save_regions_as_groups('regions')

    pb.set_bcs(ebcs=Conditions([fix_u, shift_u]))

    pb.set_solver(nls)

    status = IndexedStruct()
    state = pb.solve(status=status)

    print('Nonlinear solver status:\n', nls_status)
    print('Stationary solver status:\n', status)

    pb.save_state('linear_elasticity.vtk', state)

    if options.show:
        view = Viewer('linear_elasticity.vtk')
        view(vector_mode='warp_norm', rel_scaling=2,
             is_scalar_bar=True, is_wireframe=True)

    def _test_single_term(self, term_cls, domain, rname):
        from sfepy.terms import Term
        from sfepy.terms.terms import get_arg_kinds

        ok = True

        term_call = term_cls.name + '(%s)'

        arg_shapes_list = term_cls.arg_shapes
        if not isinstance(arg_shapes_list, list):
            arg_shapes_list = [arg_shapes_list]

        if term_cls.integration != 'custom':
            integral = self.integral

        else:
            integral = self.custom_integral

        poly_space_base = getattr(term_cls, 'poly_space_base', 'lagrange')

        prev_shapes = {}
        for _arg_shapes in arg_shapes_list:
            # Unset shapes are taken from the previous iteration.
            arg_shapes = copy(prev_shapes)
            arg_shapes.update(_arg_shapes)
            prev_shapes = arg_shapes

            self.report('arg_shapes:', arg_shapes)
            arg_types = term_cls.arg_types
            if not isinstance(arg_types[0], tuple):
                arg_types = (arg_types,)

            for iat, ats in enumerate(arg_types):
                self.report('arg_types:', ats)

                arg_kinds = get_arg_kinds(ats)
                modes = getattr(term_cls, 'modes', None)
                mode = modes[iat] if modes is not None else None

                if 'dw_s_dot_grad_i_s' in term_cls.name:
                    material_value = 0.0

                else:
                    material_value = 1.0
                aux = make_term_args(arg_shapes, arg_kinds, ats, mode, domain,
                                     material_value=material_value,
                                     poly_space_base=poly_space_base)
                args, str_args, materials, variables = aux

                self.report('args:', str_args)

                name = term_call % (', '.join(str_args))
                term = Term.new(name, integral, domain.regions[rname], **args)
                term.setup()

                call_mode = 'weak' if term.names.virtual else 'eval'
                self.report('call mode:', call_mode)

                out = term.evaluate(mode=call_mode, ret_status=True)

                if call_mode == 'eval':
                    vals, status = out
                    vals = nm.array(vals)

                else:
                    vals, iels, status = out

                if isinstance(vals, tuple):
                    # Dynamic connectivity terms.
                    vals = vals[0]

                _ok = nm.isfinite(vals).all()
                ok = ok and _ok
                self.report('values shape: %s' % (vals.shape,))
                if not _ok:
                    self.report('values are not finite!')
                    self.report(vals)

                _ok = status == 0
                if not _ok:
                    self.report('status is %d!' % status)

                ok = ok and _ok

                if term.names.virtual:
                    # Test differentiation w.r.t. state variables in the weak
                    # mode.
                    svars = term.get_state_variables(unknown_only=True)
                    for svar in svars:
                        vals, iels, status = term.evaluate(mode=call_mode,
                                                           diff_var=svar.name,
                                                           ret_status=True)
                        if isinstance(vals, tuple):
                            # Dynamic connectivity terms.
                            vals = vals[0]

                        _ok = status == 0
                        ok = ok and _ok
                        self.report('diff: %s' % svar.name)
                        if not _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 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

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(
        "-n",
        "--max-order",
        metavar="order",
        type=int,
        action="store",
        dest="max_order",
        default=10,
        help=help["max_order"],
    )
    parser.add_option(
        "-m",
        "--matrix",
        metavar="type",
        action="store",
        dest="matrix_type",
        default="laplace",
        help=help["matrix_type"],
    )
    parser.add_option(
        "-g", "--geometry", metavar="name", action="store", dest="geometry", default="2_4", help=help["geometry"]
    )
    options, args = parser.parse_args()

    dim, n_ep = int(options.geometry[0]), int(options.geometry[2])
    output("reference element geometry:")
    output("  dimension: %d, vertices: %d" % (dim, n_ep))

    n_c = {"laplace": 1, "elasticity": dim}[options.matrix_type]

    output("matrix type:", options.matrix_type)
    output("number of variable components:", n_c)

    output("polynomial space:", options.basis)

    output("max. order:", options.max_order)

    mesh = Mesh.from_file(data_dir + "/meshes/elements/%s_1.mesh" % options.geometry)
    domain = Domain("domain", mesh)
    omega = domain.create_region("Omega", "all")

    orders = nm.arange(1, options.max_order + 1, dtype=nm.int)
    conds = []

    order_fix = 0 if options.geometry in ["2_4", "3_8"] else 1

    for order in orders:
        output("order:", order, "...")

        field = Field.from_args(
            "fu", nm.float64, n_c, omega, approx_order=order, space="H1", poly_space_base=options.basis
        )

        to = field.approx_order
        quad_order = 2 * (max(to - order_fix, 0))
        output("quadrature order:", quad_order)

        integral = Integral("i", order=quad_order)
        qp, _ = integral.get_qp(options.geometry)
        output("number of quadrature points:", qp.shape[0])

        u = FieldVariable("u", "unknown", field, n_c)
        v = FieldVariable("v", "test", field, n_c, primary_var_name="u")

        m = Material("m", lam=1.0, mu=1.0)

        if options.matrix_type == "laplace":
            term = Term.new("dw_laplace(m.mu, v, u)", integral, omega, m=m, v=v, u=u)
            n_zero = 1

        else:
            assert_(options.matrix_type == "elasticity")
            term = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v, u)", integral, omega, m=m, v=v, u=u)
            n_zero = (dim + 1) * dim / 2

        term.setup()

        output("assembling...")
        tt = time.clock()
        mtx, iels = term.evaluate(mode="weak", diff_var="u")
        output("...done in %.2f s" % (time.clock() - tt))
        mtx = mtx[0][0, 0]

        try:
            assert_(nm.max(nm.abs(mtx - mtx.T)) < 1e-10)

        except:
            from sfepy.base.base import debug

            debug()

        output("matrix shape:", mtx.shape)

        eigs = eig(mtx, method="eig.sgscipy", eigenvectors=False)
        eigs.sort()
#.........这里部分代码省略.........

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()

    mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh')
    domain = Domain('domain', mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:,0]
    eps = 1e-8 * (max_x - min_x)
    omega = domain.create_region('Omega', 'all')
    gamma1 = domain.create_region('Gamma1',
                                  'nodes in x < %.10f' % (min_x + eps))
    gamma2 = domain.create_region('Gamma2',
                                  'nodes in x > %.10f' % (max_x - eps))

    field = Field('fu', nm.float64, 'vector', omega,
                  space='H1', poly_space_base='lagrange', approx_order=2)

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

    m = Material('m', lam=1.0, mu=1.0)
    f = Material('f', val=[[0.02], [0.01]])

    integral = Integral('i', order=3)

    t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u)',
         integral, omega, m=m, 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})

    bc_fun = Function('shift_u_fun', shift_u_fun, extra_args={'shift' : 0.01})
    shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = ProblemDefinition('elasticity', equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups('regions')

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print nls_status

    pb.save_state('linear_elasticity.vtk', vec)

    if options.show:
        view = Viewer('linear_elasticity.vtk')
        view(vector_mode='warp_norm', rel_scaling=2,
             is_scalar_bar=True, is_wireframe=True)

def main():
    from sfepy import data_dir

    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('--diffusivity', metavar='float', type=float,
                      action='store', dest='diffusivity',
                      default=1e-5, help=helps['diffusivity'])
    parser.add_option('--ic-max', metavar='float', type=float,
                      action='store', dest='ic_max',
                      default=2.0, help=helps['ic_max'])
    parser.add_option('--order', metavar='int', type=int,
                      action='store', dest='order',
                      default=2, help=helps['order'])
    parser.add_option('-r', '--refine', metavar='int', type=int,
                      action='store', dest='refine',
                      default=0, help=helps['refine'])
    parser.add_option('-p', '--probe',
                      action="store_true", dest='probe',
                      default=False, help=helps['probe'])
    parser.add_option('-s', '--show',
                      action="store_true", dest='show',
                      default=False, help=helps['show'])
    options, args = parser.parse_args()

    assert_((0 < options.order),
            'temperature approximation order must be at least 1!')

    output('using values:')
    output('  diffusivity:', options.diffusivity)
    output('  max. IC value:', options.ic_max)
    output('uniform mesh refinement level:', options.refine)

    mesh = Mesh.from_file(data_dir + '/meshes/3d/cylinder.mesh')
    domain = FEDomain('domain', mesh)

    if options.refine > 0:
        for ii in range(options.refine):
            output('refine %d...' % ii)
            domain = domain.refine()
            output('... %d nodes %d elements'
                   % (domain.shape.n_nod, domain.shape.n_el))

    omega = domain.create_region('Omega', 'all')
    left = domain.create_region('Left',
                                'vertices in x < 0.00001', 'facet')
    right = domain.create_region('Right',