Python Problem.time_update方法代码示例

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
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

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
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

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
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

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
    def test_solving(self):
        from sfepy.base.base import IndexedStruct
        from sfepy.discrete import (FieldVariable, Material, Problem, Function,
                                    Equation, Equations, Integral)
        from sfepy.discrete.conditions import Conditions, EssentialBC
        from sfepy.terms import Term
        from sfepy.solvers.ls import ScipyDirect
        from sfepy.solvers.nls import Newton
        from sfepy.mechanics.matcoefs import stiffness_from_lame

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

        m = Material('m', D=stiffness_from_lame(self.dim, 1.0, 1.0))
        f = Material('f', val=[[0.02], [0.01]])

        bc_fun = Function('fix_u_fun', fix_u_fun,
                          extra_args={'extra_arg' : 'hello'})

        fix_u = EssentialBC('fix_u', self.gamma1, {'u.all' : bc_fun})
        shift_u = EssentialBC('shift_u', self.gamma2, {'u.0' : 0.1})

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

        t1 = Term.new('dw_lin_elastic(m.D, v, u)',
                      integral, self.omega, m=m, v=v, u=u)

        t2 = Term.new('dw_volume_lvf(f.val, v)', integral, self.omega, f=f, v=v)

        eq = Equation('balance', t1 + t2)
        eqs = Equations([eq])

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

        state = pb.solve()

        name = op.join(self.options.out_dir, 'test_high_level_solving.vtk')
        pb.save_state(name, state)

        ok = nls_status.condition == 0
        if not ok:
            self.report('solver did not converge!')

        _ok = state.has_ebc()
        if not _ok:
            self.report('EBCs violated!')

        ok = ok and _ok

        return ok

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
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)

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
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, nls=nls, ls=ls)
    pb.time_update(ebcs=Conditions([fix_u]))

    state = pb.solve()

    return pb, state, u, gamma2

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
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

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
    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

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
def make_h1_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:`H^1` 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)
    lhs1 = Term.new('dw_volume_dot(v, %s)' % un, integral,
                    target.field.region, v=v, **{un : target})
    lhs2 = Term.new('dw_laplace(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, 'val', **kwargs)
            gval = eval_data(ts, coors, mode, 'grad', **kwargs)
            return {'val' : val, 'gval' : gval}

    m = Material('m', function=_eval_data)
    rhs1 = Term.new('dw_volume_lvf(m.val, v)', integral, target.field.region,
                    m=m, v=v)
    rhs2 = Term.new('dw_diffusion_r(m.gval, v)', integral, target.field.region,
                    m=m, v=v)

    eq = Equation('projection', lhs1 + lhs2 - rhs1 - rhs2)
    eqs = Equations([eq])

    ls = ScipyDirect({})

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

    pb = Problem('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('H1 projection: solver did not converge!')

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]

#.........这里部分代码省略.........
                        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)

    pb.time_update(ebcs=Conditions([fix1, fix2]))
    state = pb.solve()

    if subs is not None:
        field.restore_dofs()

    filename = os.path.join(options.output_dir, 'hanging.vtk')
    ensure_path(filename)

    pb.save_state(filename, state)
    if options.order > 1:
        pb.save_state(filename, state, linearization=Struct(kind='adaptive',
                                                            min_level=0,
                                                            max_level=8,
                                                            eps=1e-3))

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]

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



   
   
   # conf = Struct(method='bcgsl', precond='jacobi', sub_precond=None,
   #                 i_max=10000, eps_a=1e-50, eps_r=1e-8, eps_d=1e4,
#              verbose=True)
            
                 
                 
#ls = PETScKrylovSolver(conf)



#ls = ScipyIterative({'method':'bicgstab','i_max':100,'eps_r':1e-10})


    nls_status = IndexedStruct()
    nls = Newton({'i_max':1,'eps_a':1e-10}, lin_solver=ls, status=nls_status)

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


    dd=pb.get_materials()['Mat']
    dd.set_function(get_mat1)
    
    
    pb.save_regions_as_groups('regions')

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

    pb.time_update(ebcs=Conditions(bc))
    pb.save_regions_as_groups('regions')

#ls = ScipyIterative({'method':'bicgstab','i_max':100,'eps_r':1e-10})


#   A = pb.mtx_a
#   M = spilu(A,fill_factor = 1)
    
    #conf = Struct(solvers ='ScipyIterative',method='bcgsl', sub_precond=None,
# i_max=1000, eps_r=1e-8)
        
#pb.set_conf_solvers(conf)
    vec = pb.solve()
    print nls_status
    file.write('TIME TO SOLVE\n')
    file.write(str(nls.status.time_stats['solve'])+'\n')
    file.write('TIME TO CREATE MATRIX\n')
    file.write(str(nls.status.time_stats['matrix'])+'\n')
    #out = post_process(out, pb, state, extend=False)
    ev = pb.evaluate
    out = vec.create_output_dict()
    strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
    stress = ev('ev_cauchy_stress.3.Omega(Mat.D, u)', mode='el_avg',
                copy_materials=False)

    out['cauchy_strain'] = Struct(name='output_data', mode='cell',
                                  data=strain, dofs=None)
    out['cauchy_stress'] = Struct(name='output_data', mode='cell',
                                  data=stress, dofs=None)


    # Postprocess the solution.

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]

#.........这里部分代码省略.........
    ax = 'xyz'[:dim][options.axis]

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

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

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

    mtx_d = stiffness_from_youngpoisson(dim, options.young, options.poisson)

    m = Material('m', D=mtx_d, rho=options.density)

    integral = Integral('i', order=2*options.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)

    if options.bc_kind == 'free':
        pb.time_update()
        n_rbm = dim * (dim + 1) / 2

    elif options.bc_kind == 'cantilever':
        fixed = EssentialBC('Fixed', bottom, {'u.all' : 0.0})
        pb.time_update(ebcs=Conditions([fixed]))
        n_rbm = 0

    elif options.bc_kind == 'fixed':
        fixed = EssentialBC('Fixed', bottom_top, {'u.all' : 0.0})
        pb.time_update(ebcs=Conditions([fixed]))
        n_rbm = 0

    else:
        raise ValueError('unsupported BC kind! (%s)' % options.bc_kind)

    if options.ignore is not None:
        n_rbm = options.ignore

    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)

    try:
        eigs, svecs = eig_solver(mtx_k, mtx_m, options.n_eigs + n_rbm,
                                 eigenvectors=True)

    except sla.ArpackNoConvergence as ee:
        eigs = ee.eigenvalues

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

    mesh = omega_gi.domain.mesh

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    return pb

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]
def create_local_problem(omega_gi, order):
    """
    Local problem definition using a domain corresponding to the global region
    `omega_gi`.
    """
    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)

    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)

    field_i = Field.from_args('fu', nm.float64, 1, omega_i,
                              approx_order=order)

    output('number of local field DOFs:', field_i.n_nod)

    u_i = FieldVariable('u_i', 'unknown', field_i)
    v_i = FieldVariable('v_i', 'test', field_i, primary_var_name='u_i')

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

    mat = Material('m', lam=10, mu=5)
    t1 = Term.new('dw_laplace(m.lam, v_i, u_i)',
                  integral, omega_i, m=mat, v_i=v_i, u_i=u_i)

    def _get_load(coors):
        val = nm.ones_like(coors[:, 0])
        for coor in coors.T:
            val *= nm.sin(4 * nm.pi * coor)
        return val

    def get_load(ts, coors, mode=None, **kwargs):
        if mode == 'qp':
            return {'val' : _get_load(coors).reshape(coors.shape[0], 1, 1)}

    load = Material('load', function=Function('get_load', get_load))

    t2 = Term.new('dw_volume_lvf(load.val, v_i)',
                  integral, omega_i, load=load, v_i=v_i)

    eq = Equation('balance', t1 - 100 * t2)
    eqs = Equations([eq])

    ebc1 = EssentialBC('ebc1', gamma1_i, {'u_i.all' : 0.0})
    ebc2 = EssentialBC('ebc2', gamma2_i, {'u_i.all' : 0.1})

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

    return pb

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import time_update [as 别名]

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

    t1 = Term.new('dw_diffusion(m.diffusivity, s, T)',
                  integral, omega, m=m, s=s, T=T)
    t2 = Term.new('dw_volume_dot(s, dT/dt)',
                  integral, omega, s=s, T=T)
    eq = Equation('balance', t1 + t2)
    eqs = Equations([eq])

    # Boundary conditions.
    ebc1 = EssentialBC('T1', left, {'T.0' : 2.0})
    ebc2 = EssentialBC('T2', right, {'T.0' : -2.0})

    # Initial conditions.
    def get_ic(coors, ic):
        x, y, z = coors.T
        return 2 - 40.0 * x + options.ic_max * nm.sin(4 * nm.pi * x / 0.1)
    ic_fun = Function('ic_fun', get_ic)
    ic = InitialCondition('ic', omega, {'T.0' : ic_fun})

    pb = Problem('heat', equations=eqs)
    pb.set_bcs(ebcs=Conditions([ebc1, ebc2]))
    pb.set_ics(Conditions([ic]))

    state0 = pb.get_initial_state()
    init_fun, prestep_fun, _poststep_fun = pb.get_tss_functions(state0)

    ls = ScipyDirect({})
    nls_status = IndexedStruct()
    nls = Newton({'is_linear' : True}, lin_solver=ls, status=nls_status)
    tss = SimpleTimeSteppingSolver({'t0' : 0.0, 't1' : 100.0, 'n_step' : 11},
                                   nls=nls, context=pb, verbose=True)
    pb.set_solver(tss)

    if options.probe:
        # Prepare probe data.
        probes, labels = gen_probes(pb)

        ev = pb.evaluate
        order = 2 * (options.order - 1)

        gfield = Field.from_args('gu', nm.float64, 'vector', omega,
                                approx_order=options.order - 1)
        dvel = FieldVariable('dvel', 'parameter', gfield,
                             primary_var_name='(set-to-None)')
        cfield = Field.from_args('gu', nm.float64, 'scalar', omega,
                                approx_order=options.order - 1)
        component = FieldVariable('component', 'parameter', cfield,
                                  primary_var_name='(set-to-None)')

        nls_options = {'eps_a' : 1e-16, 'i_max' : 1}

        suffix = tss.ts.suffix
        def poststep_fun(ts, vec):
            _poststep_fun(ts, vec)

            # Probe the solution.
            dvel_qp = ev('ev_diffusion_velocity.%d.Omega(m.diffusivity, T)'
                         % order, copy_materials=False, mode='qp')
            project_by_component(dvel, dvel_qp, component, order,
                                 nls_options=nls_options)

            all_results = []
            for ii, probe in enumerate(probes):
                fig, results = probe_results(ii, T, dvel, probe, labels[ii])

                all_results.append(results)

            plt.tight_layout()
            fig.savefig('time_poisson_interactive_probe_%s.png'
                        % (suffix % ts.step), bbox_inches='tight')

            for ii, results in enumerate(all_results):
                output('probe %d (%s):' % (ii, probes[ii].name))
                output.level += 2
                for key, res in ordered_iteritems(results):
                    output(key + ':')
                    val = res[1]
                    output('  min: %+.2e, mean: %+.2e, max: %+.2e'
                           % (val.min(), val.mean(), val.max()))
                output.level -= 2

    else:
        poststep_fun = _poststep_fun

    pb.time_update(tss.ts)
    state0.apply_ebc()

    # This is required if {'is_linear' : True} is passed to Newton.
    mtx = prepare_matrix(pb, state0)
    pb.try_presolve(mtx)

    tss_status = IndexedStruct()
    tss(state0.get_vec(pb.active_only),
        init_fun=init_fun, prestep_fun=prestep_fun, poststep_fun=poststep_fun,
        status=tss_status)

    output(tss_status)

    if options.show:
        plt.show()