Python Problem.evaluate方法代码示例

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

#.........这里部分代码省略.........
    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
        svecs = ee.eigenvectors
        output('only %d eigenvalues converged!' % len(eigs))

    output('%d eigenvalues converged (%d ignored as rigid body modes)' %
           (len(eigs), n_rbm))

    eigs = eigs[n_rbm:]
    svecs = svecs[:, n_rbm:]

    omegas = nm.sqrt(eigs)
    freqs = omegas / (2 * nm.pi)

    output('number |         eigenvalue |  angular frequency '
           '|          frequency')
    for ii, eig in enumerate(eigs):
        output('%6d | %17.12e | %17.12e | %17.12e'
               % (ii + 1, eig, omegas[ii], freqs[ii]))

    # Make full eigenvectors (add DOFs fixed by boundary conditions).
    variables = pb.get_variables()

    vecs = nm.empty((variables.di.ptr[-1], svecs.shape[1]),

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import evaluate [as 别名]
    def _solve(self, property_array):
        """
        Solve the Sfepy problem for one sample.

        Args:
          property_array: array of shape (n_x, n_y, 2) where the last
          index is for Lame's parameter and shear modulus,
          respectively.

        Returns:
          the strain field of shape (n_x, n_y, 2) where the last
          index represents the x and y displacements

        """
        shape = property_array.shape[:-1]
        mesh = self._get_mesh(shape)
        domain = Domain('domain', mesh)

        region_all = domain.create_region('region_all', 'all')

        field = Field.from_args('fu', np.float64, 'vector', region_all, # pylint: disable=no-member
                                approx_order=2)

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

        m = self._get_material(property_array, domain)

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

        t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u)',
                      integral, region_all, m=m, v=v, u=u)
        eq = Equation('balance_of_forces', t1)
        eqs = Equations([eq])

        epbcs, functions = self._get_periodicBCs(domain)
        ebcs = self._get_displacementBCs(domain)
        lcbcs = self._get_linear_combinationBCs(domain)

        ls = ScipyDirect({})

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

        pb.time_update(
            ebcs=ebcs, epbcs=epbcs, lcbcs=lcbcs, functions=functions)

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

        try:
            pb.set_solvers_instances(ls, nls)
        except AttributeError:
            pb.set_solver(nls)

        vec = pb.solve()

        u = vec.create_output_dict()['u'].data
        u_reshape = np.reshape(u, (tuple(x + 1 for x in shape) + u.shape[-1:]))

        dims = domain.get_mesh_bounding_box().shape[1]
        strain = np.squeeze(
            pb.evaluate(
                'ev_cauchy_strain.{dim}.region_all(u)'.format(
                    dim=dims),
                mode='el_avg',
                copy_materials=False))
        strain_reshape = np.reshape(strain, (shape + strain.shape[-1:]))

        stress = np.squeeze(
            pb.evaluate(
                'ev_cauchy_stress.{dim}.region_all(m.D, u)'.format(
                    dim=dims),
                mode='el_avg',
                copy_materials=False))
        stress_reshape = np.reshape(stress, (shape + stress.shape[-1:]))

        return strain_reshape, u_reshape, stress_reshape