本文整理汇总了Python中sfepy.solvers.eig函数的典型用法代码示例。如果您正苦于以下问题：Python eig函数的具体用法？Python eig怎么用？Python eig使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了eig函数的11个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的Python代码示例。
示例1: solve_pressure_eigenproblem
def solve_pressure_eigenproblem(self, mtx, eig_problem=None,
n_eigs=0, check=False):
"""G = B*AI*BT or B*AI*BT+D"""
def get_slice(n_eigs, nn):
if n_eigs > 0:
ii = slice(0, n_eigs)
elif n_eigs < 0:
ii = slice(nn + n_eigs, nn)
else:
ii = slice(0, 0)
return ii
eig_problem = get_default(eig_problem, self.eig_problem)
n_eigs = get_default(n_eigs, self.n_eigs)
check = get_default(check, self.check)
mtx_c, mtx_b, action_aibt = mtx['C'], mtx['B'], mtx['action_aibt']
mtx_g = mtx_b * action_aibt.to_array() # mtx_b must be sparse!
if eig_problem == 'B*AI*BT+D':
mtx_g += mtx['D'].toarray()
mtx['G'] = mtx_g
output(mtx_c.shape, mtx_g.shape)
eigs, mtx_q = eig(mtx_c.toarray(), mtx_g, method='eig.sgscipy')
if check:
ee = nm.diag(sc.dot(mtx_q.T * mtx_c, mtx_q)).squeeze()
oo = nm.diag(sc.dot(sc.dot(mtx_q.T, mtx_g), mtx_q)).squeeze()
try:
assert_(nm.allclose(ee, eigs))
assert_(nm.allclose(oo, nm.ones_like(eigs)))
except ValueError:
debug()
nn = mtx_c.shape[0]
if isinstance(n_eigs, tuple):
output('required number of eigenvalues: (%d, %d)' % n_eigs)
if sum(n_eigs) < nn:
ii0 = get_slice(n_eigs[0], nn)
ii1 = get_slice(-n_eigs[1], nn)
eigs = nm.concatenate((eigs[ii0], eigs[ii1]))
mtx_q = nm.concatenate((mtx_q[:,ii0], mtx_q[:,ii1]), 1)
else:
output('required number of eigenvalues: %d' % n_eigs)
if (n_eigs != 0) and (abs(n_eigs) < nn):
ii = get_slice(n_eigs, nn)
eigs = eigs[ii]
mtx_q = mtx_q[:,ii]
## from sfepy.base.plotutils import pylab, iplot
## pylab.semilogy(eigs)
## pylab.figure(2)
## iplot(eigs)
## pylab.show()
## debug()
out = Struct(eigs=eigs, mtx_q=mtx_q)
return out
开发者ID:mikegraham，项目名称:sfepy，代码行数:60，代码来源:coefs_base.py
示例2: __call__
def __call__(self, problem=None, data=None):
problem = get_default(problem, self.problem)
opts = self.app_options
problem.set_equations(self.equations)
problem.select_bcs(ebc_names=self.ebcs, epbc_names=self.epbcs,
lcbc_names=self.get('lcbcs', []))
problem.update_materials(problem.ts)
self.init_solvers(problem)
mtx_a, mtx_m, data = self.prepare_matrices(problem)
output('computing resonance frequencies...')
tt = [0]
if isinstance(mtx_a, sc.sparse.spmatrix):
mtx_a = mtx_a.toarray()
if isinstance(mtx_m, sc.sparse.spmatrix):
mtx_m = mtx_m.toarray()
eigs, mtx_s_phi = eig(mtx_a, mtx_m, return_time=tt,
method=opts.eigensolver)
eigs[eigs<0.0] = 0.0
output('...done in %.2f s' % tt[0])
output('original eigenfrequencies:')
output(eigs)
opts = self.app_options
epsilon2 = opts.scale_epsilon * opts.scale_epsilon
eigs_rescaled = (opts.elasticity_contrast / epsilon2) * eigs
output('rescaled eigenfrequencies:')
output(eigs_rescaled)
output('number of eigenfrequencies: %d' % eigs.shape[0])
try:
assert_(nm.isfinite(eigs).all())
except ValueError:
from sfepy.base.base import debug; debug()
mtx_phi, eig_vectors = self.post_process(eigs, mtx_s_phi, data,
problem)
self.save(eigs, mtx_phi, problem)
evp = Struct(name='evp', eigs=eigs, eigs_rescaled=eigs_rescaled,
eig_vectors=eig_vectors)
return evp
开发者ID:Gkdnz，项目名称:sfepy，代码行数:48，代码来源:coefs_phononic.py
示例3: compute_phase_velocity
def compute_phase_velocity( self ):
from sfepy.homogenization.phono import compute_density_volume_info
opts = self.app_options
dim = self.problem.domain.mesh.dim
christoffel = self.compute_cat()
self.problem.update_materials()
dv_info = compute_density_volume_info( self.problem, opts.volume,
opts.region_to_material )
output( 'average density:', dv_info.average_density )
eye = nm.eye( dim, dim, dtype = nm.float64 )
mtx_mass = eye * dv_info.average_density
meigs, mvecs = eig( mtx_mass, mtx_b = christoffel,
eigenvectors = True, method = opts.eigensolver )
phase_velocity = 1.0 / nm.sqrt( meigs )
return phase_velocity
开发者ID:olivierverdier，项目名称:sfepy，代码行数:20，代码来源:eigen.py
示例4: main
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()
#.........这里部分代码省略.........
开发者ID:mikegraham，项目名称:sfepy，代码行数:101，代码来源:plot_condition_numbers.py
示例5: main
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 = FEDomain('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)
v = FieldVariable('v', 'test', field, 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()
# Zero 'true' zeros.
eigs[:n_zero] = 0.0
ii = nm.where(eigs < 0.0)[0]
if len(ii):
output('matrix is not positive semi-definite!')
ii = nm.where(eigs[n_zero:] < 1e-12)[0]
if len(ii):
output('matrix has more than %d zero eigenvalues!' % n_zero)
#.........这里部分代码省略.........
开发者ID:cheon7886，项目名称:sfepy，代码行数:101，代码来源:plot_condition_numbers.py
示例6: trace_full_callback
def trace_full_callback(f):
meigs, mvecs = eig((f**2) * mass(f), mtx_b=mtx_b,
eigenvectors=True, method=method)
return meigs, mvecs
开发者ID:Gkdnz，项目名称:sfepy，代码行数:5，代码来源:coefs_phononic.py
示例7: trace_callback
def trace_callback(f):
meigs = eig(mass(f), eigenvectors=False, method=method)
return meigs,
开发者ID:Gkdnz，项目名称:sfepy，代码行数:3，代码来源:coefs_phononic.py
示例8: find_zero_full_callback
def find_zero_full_callback(f):
meigs = eig((f**2) * mass(f), mtx_b=mtx_b,
eigenvectors=False, method=method)
return meigs
开发者ID:Gkdnz，项目名称:sfepy，代码行数:4，代码来源:coefs_phononic.py
示例9: find_zero_callback
def find_zero_callback(f):
meigs = eig(mass(f), eigenvectors=False, method=method)
return meigs
开发者ID:Gkdnz，项目名称:sfepy，代码行数:3，代码来源:coefs_phononic.py
示例10: solve_eigen_problem_1
def solve_eigen_problem_1( self ):
from sfepy.fem import Mesh
options = self.options
opts = self.app_options
pb = self.problem
dim = pb.domain.mesh.dim
pb.set_equations( pb.conf.equations )
pb.time_update()
output( 'assembling lhs...' )
tt = time.clock()
mtx_a = pb.evaluate(pb.conf.equations['lhs'], mode='weak',
auto_init=True, dw_mode='matrix')
output( '...done in %.2f s' % (time.clock() - tt) )
output( 'assembling rhs...' )
tt = time.clock()
mtx_b = pb.evaluate(pb.conf.equations['rhs'], mode='weak',
dw_mode='matrix')
output( '...done in %.2f s' % (time.clock() - tt) )
n_eigs = get_default( opts.n_eigs, mtx_a.shape[0] )
## mtx_a.save( 'a.txt', format='%d %d %.12f\n' )
## mtx_b.save( 'b.txt', format='%d %d %.12f\n' )
output( 'computing resonance frequencies...' )
eig = Solver.any_from_conf( pb.get_solver_conf( opts.eigen_solver ) )
eigs, mtx_s_phi = eig( mtx_a, mtx_b, n_eigs )
output( '...done' )
bounding_box = pb.domain.mesh.get_bounding_box()
# this assumes a box (3D), or a square (2D):
a = bounding_box[1][0] - bounding_box[0][0]
E_exact = None
if options.hydrogen or options.boron:
if options.hydrogen:
Z = 1
elif options.boron:
Z = 5
if dim == 2:
E_exact = [-float(Z)**2/2/(n-0.5)**2/4
for n in [1]+[2]*3+[3]*5 + [4]*8 + [5]*15]
elif dim == 3:
E_exact = [-float(Z)**2/2/n**2 for n in [1]+[2]*2**2+[3]*3**2 ]
if options.well:
if dim == 2:
E_exact = [pi**2/(2*a**2)*x
for x in [2, 5, 5, 8, 10, 10, 13, 13,
17, 17, 18, 20, 20 ] ]
elif dim == 3:
E_exact = [pi**2/(2*a**2)*x
for x in [3, 6, 6, 6, 9, 9, 9, 11, 11,
11, 12, 14, 14, 14, 14, 14,
14, 17, 17, 17] ]
if options.oscillator:
if dim == 2:
E_exact = [1] + [2]*2 + [3]*3 + [4]*4 + [5]*5 + [6]*6
elif dim == 3:
E_exact = [float(1)/2+x for x in [1]+[2]*3+[3]*6+[4]*10 ]
if E_exact is not None:
output("a=%f" % a)
output("Energies:")
output("n exact FEM error")
for i, e in enumerate(eigs):
from numpy import NaN
if i < len(E_exact):
exact = E_exact[i]
err = 100*abs((exact - e)/exact)
else:
exact = NaN
err = NaN
output("%d: %.8f %.8f %5.2f%%" % (i, exact, e, err))
else:
output(eigs)
## import sfepy.base.plotutils as plu
## plu.spy( mtx_b, eps = 1e-12 )
## plu.plt.show()
## pause()
mtx_phi = self.make_full( mtx_s_phi )
self.save_results( eigs, mtx_phi )
return Struct( pb = pb, eigs = eigs, mtx_phi = mtx_phi )
开发者ID:olivierverdier，项目名称:sfepy，代码行数:87，代码来源:schroedinger.py
示例11: solve_eigen_problem
def solve_eigen_problem( self, ofn_trunk = None, post_process_hook = None ):
if self.cached_evp is not None:
return self.cached_evp
problem = self.problem
ofn_trunk = get_default( ofn_trunk, problem.ofn_trunk,
'output file name trunk missing!' )
post_process_hook = get_default( post_process_hook,
self.post_process_hook )
conf = self.conf
eig_problem = self.app_options.eig_problem
if eig_problem in ['simple', 'simple_liquid']:
problem.set_equations( conf.equations )
problem.time_update()
mtx_a = problem.evaluate(conf.equations['lhs'], mode='weak',
auto_init=True, dw_mode='matrix')
mtx_m = problem.evaluate(conf.equations['rhs'], mode='weak',
dw_mode='matrix')
elif eig_problem == 'schur':
# A = K + B^T D^{-1} B.
mtx = assemble_by_blocks( conf.equations, self.problem,
ebcs = conf.ebcs,
epbcs = conf.epbcs )
problem.set_equations( conf.equations )
problem.time_update()
ls = Solver.any_from_conf( problem.ls_conf,
presolve = True, mtx = mtx['D'] )
mtx_b, mtx_m = mtx['B'], mtx['M']
mtx_dib = nm.empty( mtx_b.shape, dtype = mtx_b.dtype )
for ic in xrange( mtx_b.shape[1] ):
mtx_dib[:,ic] = ls( mtx_b[:,ic].toarray().squeeze() )
mtx_a = mtx['K'] + mtx_b.T * mtx_dib
else:
raise NotImplementedError
## from sfepy.base.plotutils import spy, plt
## spy( mtx_b, eps = 1e-12 )
## plt.show()
## mtx_a.save( 'a.txt', format='%d %d %.12f\n' )
## mtx_b.save( 'b.txt', format='%d %d %.12f\n' )
## pause()
output( 'computing resonance frequencies...' )
tt = [0]
if isinstance( mtx_a, sc.sparse.spmatrix ):
mtx_a = mtx_a.toarray()
if isinstance( mtx_m, sc.sparse.spmatrix ):
mtx_m = mtx_m.toarray()
eigs, mtx_s_phi = eig(mtx_a, mtx_m, return_time=tt,
method=self.app_options.eigensolver)
eigs[eigs<0.0] = 0.0
output( '...done in %.2f s' % tt[0] )
output( 'original eigenfrequencies:' )
output( eigs )
opts = self.app_options
epsilon2 = opts.scale_epsilon * opts.scale_epsilon
eigs_rescaled = (opts.elasticity_contrast / epsilon2) * eigs
output( 'rescaled eigenfrequencies:' )
output( eigs_rescaled )
output( 'number of eigenfrequencies: %d' % eigs.shape[0] )
try:
assert_( nm.isfinite( eigs ).all() )
except ValueError:
debug()
# B-orthogonality check.
## print nm.dot( mtx_s_phi[:,5], nm.dot( mtx_m, mtx_s_phi[:,5] ) )
## print nm.dot( mtx_s_phi[:,5], nm.dot( mtx_m, mtx_s_phi[:,0] ) )
## debug()
n_eigs = eigs.shape[0]
variables = problem.get_variables()
mtx_phi = nm.empty( (variables.di.ptr[-1], mtx_s_phi.shape[1]),
dtype = nm.float64 )
make_full = variables.make_full_vec
if eig_problem in ['simple', 'simple_liquid']:
for ii in xrange( n_eigs ):
mtx_phi[:,ii] = make_full( mtx_s_phi[:,ii] )
eig_vectors = mtx_phi
elif eig_problem == 'schur':
# Update also eliminated variables.
schur = self.app_options.schur
primary_var = schur['primary_var']
eliminated_var = schur['eliminated_var']
#.........这里部分代码省略.........
开发者ID:olivierverdier，项目名称:sfepy，代码行数:101，代码来源:eigen.py
注：本文中的sfepy.solvers.eig函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台，相关代码片段筛选自各路编程大神贡献的开源项目，源码版权归原作者所有，传播和使用请参考对应项目的License；未经允许，请勿转载。 |
客服电话
APP下载
官方微信
请发表评论