from __future__ import absolute_import
import numpy as nm
import numpy.linalg as nla
import scipy as sc
from sfepy.base.base import output, get_default, dict_to_struct, assert_, Struct
from sfepy.base.timing import Timer
from sfepy.solvers import eig, Solver
from sfepy.linalg import norm_l2_along_axis
from sfepy.discrete.evaluate import eval_equations
from sfepy.homogenization.coefs_base import MiniAppBase, CorrMiniApp
from sfepy.homogenization.utils import coor_to_sym
import six
from six.moves import range
[docs]def compute_eigenmomenta(em_equation, var_name, problem, eig_vectors,
transform=None):
"""
Compute the eigenmomenta corresponding to given eigenvectors.
"""
n_dof, n_eigs = eig_vectors.shape
equations, variables = problem.create_evaluable(em_equation)
variables.init_state()
var = variables[var_name]
n_c = var.n_components
eigenmomenta = nm.empty((n_eigs, n_c), dtype=nm.float64)
for ii in range(n_eigs):
if transform is None:
vec_phi, is_zero = eig_vectors[:,ii], False
else:
vec_phi, is_zero = transform(eig_vectors[:,ii], (n_dof / n_c, n_c))
if is_zero:
eigenmomenta[ii, :] = 0.0
else:
variables.set_state_parts({var_name: vec_phi})
val = eval_equations(equations, variables)
eigenmomenta[ii, :] = val
return eigenmomenta
[docs]def get_ranges(freq_range, eigs):
"""
Get an eigenvalue range slice and a corresponding initial frequency range
within a given frequency range.
"""
mine, maxe = freq_range
ii = nm.where((eigs > (mine**2.)) & (eigs < (maxe**2.)))[0]
freq_range_initial = nm.sqrt(eigs[ii])
eig_range = (ii[0], ii[-1] + 1) # +1 as it is a slice.
eig_range = slice(*eig_range)
return freq_range_initial, eig_range
[docs]def cut_freq_range(freq_range, eigs, valid, freq_margins, eig_range,
fixed_freq_range, freq_eps):
"""
Cut off masked resonance frequencies. Margins are preserved, like no
resonances were cut.
Returns
-------
freq_range : array
The new range of frequencies.
freq_range_margins : array
The range of frequencies with prepended/appended margins equal to
`fixed_freq_range` if it is not None.
"""
n_eigs = eigs.shape[0]
output('masked resonance frequencies in range:')
output(nm.where(valid[eig_range] == False)[0])
if fixed_freq_range is None:
min_freq, max_freq = freq_range[0], freq_range[-1]
margins = freq_margins * (max_freq - min_freq)
prev_freq = min_freq - margins[0]
next_freq = max_freq + margins[1]
if eig_range.start > 0:
prev_freq = max(nm.sqrt(eigs[eig_range.start - 1]) + freq_eps,
prev_freq)
if eig_range.stop < n_eigs:
next_freq = min(nm.sqrt(eigs[eig_range.stop]) - freq_eps,
next_freq)
prev_freq = max(freq_eps, prev_freq)
next_freq = max(freq_eps, next_freq, prev_freq + freq_eps)
else:
prev_freq, next_freq = fixed_freq_range
freq_range = freq_range[valid[eig_range]]
freq_range_margins = nm.r_[prev_freq, freq_range, next_freq]
return freq_range, freq_range_margins
[docs]def split_chunks(indx):
"""Split index vector to chunks of consecutive numbers."""
if not len(indx): return []
delta = nm.ediff1d(indx, to_end=2)
ir = nm.where(delta > 1)[0]
chunks = []
ic0 = 0
for ic in ir:
chunk = indx[ic0:ic+1]
ic0 = ic + 1
chunks.append(chunk)
return chunks
[docs]def get_log_freqs(f0, f1, df, freq_eps, n_point_min, n_point_max):
"""
Get logging frequencies.
The frequencies get denser towards the interval boundaries.
"""
f_delta = f1 - f0
f_mid = 0.5 * (f0 + f1)
if (f1 - f0) > (2.0 * freq_eps):
num = min(n_point_max, max(n_point_min, int((f1 - f0) / df)))
a = nm.linspace(0., 1., num)
log_freqs = f0 + freq_eps \
+ 0.5 * (nm.sin((a - 0.5) * nm.pi) + 1.0) \
* (f1 - f0 - 2.0 * freq_eps)
else:
log_freqs = nm.array([f_mid - 1e-8 * f_delta,
f_mid + 1e-8 * f_delta])
return log_freqs
[docs]def detect_band_gaps(mass, freq_info, opts, gap_kind='normal', mtx_b=None):
"""
Detect band gaps given solution to eigenproblem (eigs,
eig_vectors). Only valid resonance frequencies (e.i. those for which
corresponding eigenmomenta are above a given threshold) are taken into
account.
Notes
-----
- make freq_eps relative to ]f0, f1[ size?
"""
output('eigensolver:', opts.eigensolver)
fm = freq_info.freq_range_margins
min_freq, max_freq = fm[0], fm[-1]
output('freq. range with margins: [%8.3f, %8.3f]'
% (min_freq, max_freq))
df = opts.freq_step * (max_freq - min_freq)
fz_callback = get_callback(mass.evaluate, opts.eigensolver,
mtx_b=mtx_b, mode='find_zero')
trace_callback = get_callback(mass.evaluate, opts.eigensolver,
mtx_b=mtx_b, mode='trace')
n_col = 1 + (mtx_b is not None)
logs = [[] for ii in range(n_col + 1)]
gaps = []
for ii in range(freq_info.freq_range.shape[0] + 1):
f0, f1 = fm[[ii, ii+1]]
output('interval: ]%.8f, %.8f[...' % (f0, f1))
log_freqs = get_log_freqs(f0, f1, df, opts.freq_eps, 100, 1000)
output('n_logged: %d' % log_freqs.shape[0])
log_mevp = [[] for ii in range(n_col)]
for f in log_freqs:
for ii, data in enumerate(trace_callback(f)):
log_mevp[ii].append(data)
# Get log for the first and last f in log_freqs.
lf0 = log_freqs[0]
lf1 = log_freqs[-1]
log0, log1 = log_mevp[0][0], log_mevp[0][-1]
min_eig0 = log0[0]
max_eig1 = log1[-1]
if gap_kind == 'liquid':
mevp = nm.array(log_mevp, dtype=nm.float64).squeeze()
si = nm.where(mevp[:,0] < 0.0)[0]
li = nm.where(mevp[:,-1] < 0.0)[0]
wi = nm.setdiff1d(si, li)
if si.shape[0] == 0: # No gaps.
gap = ([2, lf0, log0[0]], [2, lf0, log0[-1]])
gaps.append(gap)
elif li.shape[0] == mevp.shape[0]: # Full interval strong gap.
gap = ([1, lf1, log1[0]], [1, lf1, log1[-1]])
gaps.append(gap)
else:
subgaps = []
for chunk in split_chunks(li): # Strong gaps.
i0, i1 = chunk[0], chunk[-1]
fmin, fmax = log_freqs[i0], log_freqs[i1]
gap = ([1, fmin, mevp[i0,-1]], [1, fmax, mevp[i1,-1]])
subgaps.append(gap)
for chunk in split_chunks(wi): # Weak gaps.
i0, i1 = chunk[0], chunk[-1]
fmin, fmax = log_freqs[i0], log_freqs[i1]
gap = ([0, fmin, mevp[i0,-1]], [2, fmax, mevp[i1,-1]])
subgaps.append(gap)
gaps.append(subgaps)
else:
if min_eig0 > 0.0: # No gaps.
gap = ([2, lf0, log0[0]], [2, lf0, log0[-1]])
elif max_eig1 < 0.0: # Full interval strong gap.
gap = ([1, lf1, log1[0]], [1, lf1, log1[-1]])
else:
llog_freqs = list(log_freqs)
# Insert fmin, fmax into log.
output('finding zero of the largest eig...')
smax, fmax, vmax = find_zero(lf0, lf1, fz_callback,
opts.freq_eps, opts.zero_eps, 1)
im = nm.searchsorted(log_freqs, fmax)
llog_freqs.insert(im, fmax)
for ii, data in enumerate(trace_callback(fmax)):
log_mevp[ii].insert(im, data)
output('...done')
if smax in [0, 2]:
output('finding zero of the smallest eig...')
# having fmax instead of f0 does not work if freq_eps is
# large.
smin, fmin, vmin = find_zero(lf0, lf1, fz_callback,
opts.freq_eps, opts.zero_eps, 0)
im = nm.searchsorted(log_freqs, fmin)
# +1 due to fmax already inserted before.
llog_freqs.insert(im+1, fmin)
for ii, data in enumerate(trace_callback(fmin)):
log_mevp[ii].insert(im+1, data)
output('...done')
elif smax == 1:
smin = 1 # both are negative everywhere.
fmin, vmin = fmax, vmax
gap = ([smin, fmin, vmin], [smax, fmax, vmax])
log_freqs = nm.array(llog_freqs)
output(gap[0])
output(gap[1])
gaps.append(gap)
logs[0].append(log_freqs)
for ii, data in enumerate(log_mevp):
logs[ii+1].append(nm.array(data, dtype = nm.float64))
output('...done')
kinds = describe_gaps(gaps)
slogs = Struct(freqs=logs[0], eigs=logs[1])
if n_col == 2:
slogs.eig_vectors = logs[2]
return slogs, gaps, kinds
[docs]def get_callback(mass, solver_kind, mtx_b=None, mode='trace'):
"""
Return callback to solve band gaps or dispersion eigenproblem P.
Notes
-----
Find zero callbacks return:
eigenvalues
Trace callbacks return:
(eigenvalues,)
or
(eigenvalues, eigenvectors) (in full (dispoersion) mode)
If `mtx_b` is None, the problem P is
M w = \lambda w,
otherwise it is
omega^2 M w = \eta B w"""
def find_zero_callback(f):
meigs = eig(mass(f), eigenvectors=False, solver_kind=solver_kind)
return meigs
def find_zero_full_callback(f):
meigs = eig((f**2) * mass(f), mtx_b=mtx_b,
eigenvectors=False, solver_kind=solver_kind)
return meigs
def trace_callback(f):
meigs = eig(mass(f), eigenvectors=False, solver_kind=solver_kind)
return meigs,
def trace_full_callback(f):
meigs, mvecs = eig((f**2) * mass(f), mtx_b=mtx_b,
eigenvectors=True, solver_kind=solver_kind)
return meigs, mvecs
if mtx_b is not None:
mode += '_full'
return eval(mode + '_callback')
[docs]def find_zero(f0, f1, callback, freq_eps, zero_eps, mode):
"""
For f \in ]f0, f1[ find frequency f for which either the smallest (`mode` =
0) or the largest (`mode` = 1) eigenvalue of problem P given by `callback`
is zero.
Returns
-------
flag : 0, 1, or 2
The flag, see Notes below.
frequency : float
The found frequency.
eigenvalue : float
The eigenvalue corresponding to the found frequency.
Notes
-----
Meaning of the return value combinations:
===== ====== ========
mode flag meaning
===== ====== ========
0, 1 0 eigenvalue -> 0 for f \in ]f0, f1[
0 1 f -> f1, smallest eigenvalue < 0
0 2 f -> f0, smallest eigenvalue > 0 and -> -\infty
1 1 f -> f1, largest eigenvalue < 0 and -> +\infty
1 2 f -> f0, largest eigenvalue > 0
===== ====== ========
"""
fm, fp = f0, f1
ieig = {0 : 0, 1 : -1}[mode]
while 1:
f = 0.5 * (fm + fp)
meigs = callback(f)
val = meigs[ieig]
## print f, f0, f1, fm, fp, val
## print '%.16e' % f, '%.16e' % fm, '%.16e' % fp, '%.16e' % val
if ((abs(val) < zero_eps)
or ((fp - fm) < (abs(fm) * nm.finfo(float).eps))):
return 0, f, val
if mode == 0:
if (f - f0) < freq_eps:
return 2, f0, val
elif (f1 - f) < freq_eps:
return 1, f1, val
elif mode == 1:
if (f1 - f) < freq_eps:
return 1, f1, val
elif (f - f0) < freq_eps:
return 2, f0, val
if val > 0.0:
fp = f
else:
fm = f
[docs]def describe_gaps(gaps):
kinds = []
for ii, gap in enumerate(gaps):
if isinstance(gap, list):
subkinds = []
for gmin, gmax in gap:
if (gmin[0] == 2) and (gmax[0] == 2):
kind = ('p', 'propagation zone')
elif (gmin[0] == 1) and (gmax[0] == 1):
kind = ('is', 'inner strong band gap')
elif (gmin[0] == 0) and (gmax[0] == 2):
kind = ('iw', 'inner weak band gap')
subkinds.append(kind)
kinds.append(subkinds)
else:
gmin, gmax = gap
if (gmin[0] == 2) and (gmax[0] == 2):
kind = ('p', 'propagation zone')
elif (gmin[0] == 1) and (gmax[0] == 2):
kind = ('w', 'full weak band gap')
elif (gmin[0] == 0) and (gmax[0] == 2):
kind = ('wp', 'weak band gap + propagation zone')
elif (gmin[0] == 1) and (gmax[0] == 1):
kind = ('s', 'full strong band gap (due to end of freq.'
' range or too large thresholds)')
elif (gmin[0] == 1) and (gmax[0] == 0):
kind = ('sw', 'strong band gap + weak band gap')
elif (gmin[0] == 0) and (gmax[0] == 0):
kind = ('swp', 'strong band gap + weak band gap +'
' propagation zone')
else:
msg = 'impossible band gap combination: %d, %d' % (gmin, gmax)
raise ValueError(msg)
kinds.append(kind)
return kinds
[docs]def get_gap_ranges(freq_range, gaps, kinds):
"""
For each (potential) band gap in `gaps`, return the frequency ranges of its
parts according to `kinds`.
"""
def get_ranges(ii, f0, f1, kind, kind_desc, gmin, gmax):
if kind == 'p':
ranges = [(f0, f1)]
elif kind == 'w':
ranges = [(f0, f1)]
elif kind == 'wp':
ranges = [(f0, gmin[1]), (gmin[1], f1)]
elif kind == 's':
ranges = [(f0, f1)]
elif kind == 'sw':
ranges = [(f0, gmax[1]), (gmax[1], f1)]
elif kind == 'swp':
ranges = [(f0, gmax[1]), (gmax[1], gmin[1]), (gmin[1], f1)]
elif kind == 'is':
ranges = [(gmin[1], gmax[1])]
elif kind == 'iw':
ranges = [(gmin[1], gmax[1])]
else:
msg = 'impossible band gap combination! (%f %d)' % (gmin, gmax)
raise ValueError(msg)
return ranges
gap_ranges = []
for ii in range(len(freq_range) - 1):
f0, f1 = freq_range[[ii, ii+1]]
gap = gaps[ii]
if isinstance(gap, list):
ranges = []
for ig, (gmin, gmax) in enumerate(gap):
kind, kind_desc = kinds[ii][ig]
aux = get_ranges(ii, f0, f1, kind, kind_desc, gmin, gmax)
ranges.append(aux)
else:
gmin, gmax = gap
kind, kind_desc = kinds[ii]
ranges = get_ranges(ii, f0, f1, kind, kind_desc, gmin, gmax)
gap_ranges.append(ranges)
return gap_ranges
[docs]def compute_cat_sym_sym(coef, iw_dir):
"""
Christoffel acoustic tensor (part) of elasticity tensor dimension.
"""
dim = iw_dir.shape[0]
cat = nm.zeros((dim, dim), dtype=nm.float64)
for ii in range(dim):
for ij in range(dim):
ir = coor_to_sym(ii, ij, dim)
for ik in range(dim):
for il in range(dim):
ic = coor_to_sym(ik, il, dim)
cat[ii,ik] += coef[ir,ic] * iw_dir[ij] * iw_dir[il]
return cat
[docs]def compute_cat_dim_sym(coef, iw_dir):
"""
Christoffel acoustic tensor part of piezo-coupling tensor dimension.
"""
dim = iw_dir.shape[0]
cat = nm.zeros((dim,), dtype=nm.float64)
for ii in range(dim):
for ij in range(dim):
ir = coor_to_sym(ii, ij, dim)
for ik in range(dim):
cat[ii] += coef[ik,ir] * iw_dir[ij] * iw_dir[ik]
return cat
[docs]def compute_cat_dim_dim(coef, iw_dir):
"""
Christoffel acoustic tensor part of dielectric tensor dimension.
"""
cat = nm.dot(nm.dot(coef, iw_dir), iw_dir)
return cat
[docs]class SimpleEVP(CorrMiniApp):
"""
Simple eigenvalue problem.
"""
[docs] def process_options(self):
get = self.options.get
return Struct(eigensolver=get('eigensolver', 'eig.sgscipy'),
elasticity_contrast=get('elasticity_contrast', 1.0),
scale_epsilon=get('scale_epsilon', 1.0),
save_eig_vectors=get('save_eig_vectors', (0, 0)))
def __call__(self, problem=None, data=None):
problem = get_default(problem, self.problem)
opts = self.app_options
if self.equations is not None:
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,
solver_kind=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
[docs] def prepare_matrices(self, problem):
mtx_a = problem.evaluate(self.equations['lhs'], mode='weak',
auto_init=True, dw_mode='matrix')
mtx_m = problem.evaluate(self.equations['rhs'], mode='weak',
dw_mode='matrix')
return mtx_a, mtx_m, None
[docs] def post_process(self, eigs, mtx_s_phi, data, problem):
n_eigs = eigs.shape[0]
variables = problem.get_variables()
mtx_phi = nm.empty((variables.di.n_dof_total, mtx_s_phi.shape[1]),
dtype=nm.float64)
make_full = variables.make_full_vec
for ii in range(n_eigs):
mtx_phi[:,ii] = make_full(mtx_s_phi[:,ii])
return mtx_phi, mtx_phi
[docs] def save(self, eigs, mtx_phi, problem):
save = self.app_options.save_eig_vectors
n_eigs = eigs.shape[0]
out = {}
variables = problem.set_default_state()
for ii in range(n_eigs):
if (ii >= save[0]) and (ii < (n_eigs - save[1])): continue
variables.set_state(mtx_phi[:,ii], force=True)
aux = variables.create_output()
for name, val in six.iteritems(aux):
out[name+'%03d' % ii] = val
if self.post_process_hook is not None:
out = self.post_process_hook(out, problem, mtx_phi)
problem.domain.mesh.write(self.save_name + '.vtk', io='auto', out=out)
fd = open(self.save_name + '_eigs.txt', 'w')
eigs.tofile(fd, ' ')
fd.close()
[docs]class SchurEVP(SimpleEVP):
"""
Schur complement eigenvalue problem.
"""
[docs] def prepare_matrices(self, problem):
"""
A = K + B^T D^{-1} B
"""
equations = problem.equations
mtx = equations.eval_tangent_matrices(None, problem.mtx_a,
by_blocks=True)
ls = Solver.any_from_conf(problem.ls_conf
+ Struct(use_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 range(mtx_b.shape[1]):
mtx_dib[:,ic] = ls(mtx_b[:,ic].toarray().squeeze())
mtx_a = mtx['K'] + mtx_b.T * mtx_dib
return mtx_a, mtx_m, mtx_dib
[docs] def post_process(self, eigs, mtx_s_phi, mtx_dib, problem):
n_eigs = eigs.shape[0]
variables = problem.get_variables()
mtx_phi = nm.empty((variables.di.n_dof_total, mtx_s_phi.shape[1]),
dtype=nm.float64)
make_full = variables.make_full_vec
# Update also eliminated variables.
schur = self.app_options.schur
primary_var = schur['primary_var']
eliminated_var = schur['eliminated_var']
mtx_s_phi_schur = - sc.dot(mtx_dib, mtx_s_phi)
aux = nm.empty((variables.adi.n_dof_total,), dtype=nm.float64)
setv = variables.set_vec_part
for ii in range(n_eigs):
setv(aux, primary_var, mtx_s_phi[:,ii], reduced=True)
setv(aux, eliminated_var, mtx_s_phi_schur[:,ii],
reduced=True)
mtx_phi[:,ii] = make_full(aux)
indx = variables.get_indx(primary_var)
eig_vectors = mtx_phi[indx,:]
return mtx_phi, eig_vectors
[docs]class DensityVolumeInfo(MiniAppBase):
"""
Determine densities of regions specified in `region_to_material`, and
compute average density based on region volumes.
"""
def __call__(self, volume=None, problem=None, data=None):
problem = get_default(problem, self.problem)
vf = data[self.requires[0]]
average_density = 0.0
total_volume = 0.0
volumes = {}
densities = {}
for region_name, aux in six.iteritems(self.region_to_material):
vol = vf['volume_' + region_name]
mat_name, item_name = aux
conf = problem.conf.get_item_by_name('materials', mat_name)
density = conf.values[item_name]
output('region %s: volume %f, density %f' % (region_name,
vol, density))
volumes[region_name] = vol
densities[region_name] = density
average_density += vol * density
total_volume += vol
true_volume = self._get_volume(volume)
output('total volume:', true_volume)
average_density /= true_volume
return Struct(name='density_volume_info',
average_density=average_density,
total_volume=true_volume,
volumes=volumes,
densities=densities,
to_file_txt=self.to_file_txt)
[docs] @staticmethod
def to_file_txt(fd, float_format, dv_info):
ff = float_format + '\n'
fd.write('total volume:\n')
fd.write(ff % dv_info.total_volume)
fd.write('average density:\n')
fd.write(ff % dv_info.average_density)
for key, val in six.iteritems(dv_info.volumes):
fd.write('%s volume:\n' % key)
fd.write(ff % val)
fd.write('%s density:\n' % key)
fd.write(ff % dv_info.densities[key])
[docs]class Eigenmomenta(MiniAppBase):
"""
Eigenmomenta corresponding to eigenvectors.
Parameters
----------
var_name : str
The name of the variable used in the integral.
threshold : float
The threshold under which an eigenmomentum is considered zero.
threshold_is_relative : bool
If True, the `threshold` is relative w.r.t. max. norm of eigenmomenta.
transform : callable, optional
Optional function for transforming the eigenvectors before computing
the eigenmomenta.
Returns
-------
eigenmomenta : Struct
The resulting eigenmomenta. An eigenmomentum above threshold is marked
by the attribute 'valid' set to True.
"""
[docs] def process_options(self):
options = dict_to_struct(self.options)
get = options.get
return Struct(var_name=get('var_name', None,
'missing "var_name" in options!'),
threshold=get('threshold', 1e-4),
threshold_is_relative=get('threshold_is_relative', True),
transform=get('transform', None))
def __call__(self, volume=None, problem=None, data=None):
problem = get_default(problem, self.problem)
opts = self.app_options
evp, dv_info = [data[ii] for ii in self.requires]
output('computing eigenmomenta...')
if opts.transform is not None:
fun = problem.conf.get_function(opts.transform[0])
def wrap_transform(vec, shape):
return fun(vec, shape, *opts.eig_vector_transform[1:])
else:
wrap_transform = None
timer = Timer(start=True)
eigenmomenta = compute_eigenmomenta(self.expression, opts.var_name,
problem, evp.eig_vectors,
wrap_transform)
output('...done in %.2f s' % timer.stop())
n_eigs = evp.eigs.shape[0]
mag = norm_l2_along_axis(eigenmomenta)
if opts.threshold_is_relative:
tol = opts.threshold * mag.max()
else:
tol = opts.threshold
valid = nm.where(mag < tol, False, True)
mask = nm.where(valid == False)[0]
eigenmomenta[mask, :] = 0.0
n_zeroed = mask.shape[0]
output('%d of %d eigenmomenta zeroed (under %.2e)'\
% (n_zeroed, n_eigs, tol))
out = Struct(name='eigenmomenta', n_zeroed=n_zeroed,
eigenmomenta=eigenmomenta, valid=valid,
to_file_txt=None)
return out
[docs]class AcousticMassTensor(MiniAppBase):
"""
The acoustic mass tensor for a given frequency.
Returns
-------
self : AcousticMassTensor instance
This class instance whose `evaluate()` method computes for a given
frequency the required tensor.
Notes
-----
`eigenmomenta`, `eigs` should contain only valid resonances.
"""
to_file_txt = None
def __call__(self, volume=None, problem=None, data=None):
evp, self.dv_info, ema = [data[ii] for ii in self.requires]
self.eigs = evp.eigs[ema.valid]
self.eigenmomenta = ema.eigenmomenta[ema.valid, :]
return self
[docs] def evaluate(self, freq):
ema = self.eigenmomenta
n_c = ema.shape[1]
fmass = nm.zeros((n_c, n_c), dtype=nm.float64)
num, denom = self.get_coefs(freq)
de = 1.0 / denom
if not nm.isfinite(de).all():
raise ValueError('frequency %e too close to resonance!' % freq)
for ir in range(n_c):
for ic in range(n_c):
if ir <= ic:
val = nm.sum(num * de * (ema[:, ir] * ema[:, ic]))
fmass[ir, ic] += val
else:
fmass[ir, ic] = fmass[ic, ir]
eye = nm.eye(n_c, n_c, dtype=nm.float64)
mtx_mass = (eye * self.dv_info.average_density) \
- (fmass / self.dv_info.total_volume)
return mtx_mass
[docs] def get_coefs(self, freq):
"""
Get frequency-dependent coefficients.
"""
f2 = freq*freq
de = f2 - self.eigs
return f2, de
[docs]class AcousticMassLiquidTensor(AcousticMassTensor):
[docs] def get_coefs(self, freq):
"""
Get frequency-dependent coefficients.
"""
eigs = self.eigs
f2 = freq*freq
aux = (f2 - self.gamma * eigs)
num = f2 * aux
denom = aux*aux + f2*(self.eta*self.eta)*nm.power(eigs, 2.0)
return num, denom
[docs]class AppliedLoadTensor(AcousticMassTensor):
"""
The applied load tensor for a given frequency.
Returns
-------
self : AppliedLoadTensor instance
This class instance whose `evaluate()` method computes for a given
frequency the required tensor.
Notes
-----
`eigenmomenta`, `ueigenmomenta`, `eigs` should contain only valid
resonances.
"""
to_file_txt = None
def __call__(self, volume=None, problem=None, data=None):
evp, self.dv_info, ema, uema = [data[ii] for ii in self.requires]
self.eigs = evp.eigs[ema.valid]
self.eigenmomenta = ema.eigenmomenta[ema.valid, :]
self.ueigenmomenta = uema.eigenmomenta[uema.valid, :]
return self
[docs] def evaluate(self, freq):
ema, uema = self.eigenmomenta, self.ueigenmomenta
n_c = ema.shape[1]
fload = nm.zeros((n_c, n_c), dtype=nm.float64)
num, denom = self.get_coefs(freq)
de = 1.0 / denom
if not nm.isfinite(de).all():
raise ValueError('frequency %e too close to resonance!' % freq)
for ir in range(n_c):
for ic in range(n_c):
val = nm.sum(num * de * (ema[:, ir] * uema[:, ic]))
fload[ir, ic] += val
eye = nm.eye(n_c, n_c, dtype=nm.float64)
mtx_load = eye - (fload / self.dv_info.total_volume)
return mtx_load
[docs]class BandGaps(MiniAppBase):
"""
Band gaps detection.
Parameters
----------
eigensolver : str
The name of the eigensolver for mass matrix eigenvalues.
eig_range : (int, int)
The eigenvalues range (squared frequency) to consider.
freq_margins : (float, float)
Margins in percents of initial frequency range given by
`eig_range` by which the range is increased.
fixed_freq_range : (float, float)
The frequency range to consider. Has precedence over `eig_range`
and `freq_margins`.
freq_step : float
The frequency step for tracing, in percent of the frequency range.
freq_eps : float
The frequency difference smaller than `freq_eps` is considered zero.
zero_eps : float
The tolerance for finding zeros of mass matrix eigenvalues.
detect_fun : callable
The function for detecting the band gaps. Default is
:func:`detect_band_gaps()`.
log_save_name : str
If not None, the band gaps log is to be saved under the given name.
raw_log_save_name : str
If not None, the raw band gaps log is to be saved under the given name.
"""
[docs] def process_options(self):
get = self.options.get
freq_margins = get('freq_margins', (5, 5))
# Given per cent.
freq_margins = 0.01 * nm.array(freq_margins, dtype=nm.float64)
# Given in per cent.
freq_step = 0.01 * get('freq_step', 5)
return Struct(eigensolver=get('eigensolver', 'eig.sgscipy'),
eig_range=get('eig_range', None),
freq_margins=freq_margins,
fixed_freq_range=get('fixed_freq_range', None),
freq_step=freq_step,
freq_eps=get('freq_eps', 1e-8),
zero_eps=get('zero_eps', 1e-8),
detect_fun=get('detect_fun', detect_band_gaps),
log_save_name=get('log_save_name', None),
raw_log_save_name=get('raw_log_save_name', None))
def __call__(self, volume=None, problem=None, data=None):
problem = get_default(problem, self.problem)
opts = self.app_options
evp, ema, mass = [data[ii] for ii in self.requires[:3]]
if len(self.requires) == 4:
mtx_b = data[self.requires[3]]
else:
mtx_b = None
eigs = evp.eigs
self.fix_eig_range(eigs.shape[0])
if opts.fixed_freq_range is not None:
(freq_range_initial,
opts.eig_range) = get_ranges(opts.fixed_freq_range, eigs)
else:
opts.eig_range = slice(*opts.eig_range)
freq_range_initial = nm.sqrt(eigs[opts.eig_range])
output('initial freq. range : [%8.3f, %8.3f]'
% tuple(freq_range_initial[[0, -1]]))
aux = cut_freq_range(freq_range_initial, eigs, ema.valid,
opts.freq_margins, opts.eig_range,
opts.fixed_freq_range,
opts.freq_eps)
freq_range, freq_range_margins = aux
if len(freq_range):
output('freq. range : [%8.3f, %8.3f]'
% tuple(freq_range[[0, -1]]))
else:
# All masked.
output('freq. range : all masked!')
freq_info = Struct(name='freq_info',
freq_range_initial=freq_range_initial,
freq_range=freq_range,
freq_range_margins=freq_range_margins)
logs, gaps, kinds = opts.detect_fun(mass, freq_info, opts, mtx_b=mtx_b)
gap_ranges = get_gap_ranges(freq_range_margins, gaps, kinds)
bg = Struct(name='band_gaps', logs=logs, gaps=gaps, kinds=kinds,
gap_ranges=gap_ranges,
valid=ema.valid, eig_range=opts.eig_range,
n_eigs=eigs.shape[0], n_zeroed=ema.n_zeroed,
freq_range_initial=freq_info.freq_range_initial,
freq_range=freq_info.freq_range,
freq_range_margins=freq_info.freq_range_margins,
opts=opts, to_file_txt=self.to_file_txt,
log_save_name=opts.log_save_name,
raw_log_save_name=opts.raw_log_save_name,
save_log=self.save_log)
return bg
[docs] def fix_eig_range(self, n_eigs):
eig_range = get_default(self.app_options.eig_range, (0, n_eigs))
if eig_range[-1] < 0:
eig_range[-1] += n_eigs + 1
assert_(eig_range[0] < (eig_range[1] - 1))
assert_(eig_range[1] <= n_eigs)
self.app_options.eig_range = eig_range
[docs] @staticmethod
def to_file_txt(fd, float_format, bg):
if bg.log_save_name is not None:
fd.write(bg.log_save_name + '\n')
else:
fd.write('--\n')
[docs] @staticmethod
def save_log(filename, float_format, bg):
"""
Save band gaps, valid flags and eigenfrequencies.
"""
fd = open(filename, 'w')
freq_range = bg.freq_range_margins
fd.write('n_zeroed: %d\n' % bg.n_zeroed)
fd.write('n_eigs: %d\n' % bg.n_eigs)
n_row = len(freq_range) - 1
fd.write('n_ranges: %d\n' % n_row)
fd.write('f0 f1 flag_min f_min v_min flag_max f_max v_max'
' kind\ndesc\n')
ff = float_format
format = "%s %s %%d %s %s %%d %s %s %%s\n%%s\n" % (6 * (ff,))
for ir in range(n_row):
f0, f1 = freq_range[[ir, ir+1]]
gmin, gmax = bg.gaps[ir]
fd.write(format % ((f0, f1) + tuple(gmin) + tuple(gmax)
+ bg.kinds[ir]))
fd.write('\nname kind f_from f_to index f0 f1\n')
format = '%%s %%s %s %s %%d %s %s\n' % (4 * (ff,))
for ii, f0 in enumerate(bg.freq_range_margins[:-1]):
f1 = bg.freq_range_margins[ii + 1]
kind = bg.kinds[ii][0]
for ir, rng in enumerate(bg.gap_ranges[ii]):
fd.write(format
% (bg.name, kind[ir], rng[0], rng[1], ii, f0, f1))
n_row = len(freq_range)
fd.write('\nn_resonance: %d\n' % n_row)
fd.write('valid f\n')
freq_range = bg.freq_range_initial
valid_in_range = bg.valid[bg.eig_range]
format = "%%d %s\n" % ff
for ir in range(n_row):
fd.write(format % (valid_in_range[ir], freq_range[ir]))
fd.close()
[docs]class ChristoffelAcousticTensor(MiniAppBase):
[docs] def process_options(self):
get = self.options.get
return Struct(mode=get('mode', 'simple'),
incident_wave_dir=get('incident_wave_dir', None))
r"""
Compute Christoffel acoustic tensor (cat) given the incident wave
direction (unit vector).
Parameters
----------
mode : 'simple' or 'piezo'
The call mode.
incident_wave_dir : array
The incident wave direction vector.
Returns
-------
cat : array
The Christoffel acoustic tensor.
Notes
-----
- If mode == 'simple', only the elasticity tensor :math:`C_{ijkl}` is used
and cat := :math:`\Gamma_{ik} = C_{ijkl} n_j n_l`.
- If mode == 'piezo', also the piezo-coupling :math:`G_{ijk}` and
dielectric :math:`D_{ij}` tensors are used and cat := :math:`H_{ik} =
\Gamma_{ik} + \frac{1}{\xi} \gamma_i \gamma_j`, where :math:`\gamma_i =
G_{kij} n_j n_k`, :math:`\xi = D_{kl} n_k n_l`.
"""
def __call__(self, volume=None, problem=None, data=None):
problem = get_default(problem, self.problem)
opts = self.app_options
iw_dir = nm.array(opts.incident_wave_dir, dtype=nm.float64)
dim = problem.get_dim()
assert_(dim == iw_dir.shape[0])
iw_dir = iw_dir / nla.norm(iw_dir)
elastic = data[self.requires[0]]
cat = compute_cat_sym_sym(elastic, iw_dir)
if opts.mode =='piezo':
dielectric, coupling = [data[ii] for ii in self.requires[1:]]
xi = compute_cat_dim_dim(dielectric, iw_dir)
gamma = compute_cat_dim_sym(coupling, iw_dir)
cat += nm.outer(gamma, gamma) / xi
return cat
[docs]class PolarizationAngles(MiniAppBase):
"""
Compute polarization angles, i.e., angles between incident wave direction
and wave vectors. Vector length does not matter - eigenvectors are used
directly.
"""
[docs] def process_options(self):
get = self.options.get
return Struct(incident_wave_dir=get('incident_wave_dir', None))
def __call__(self, volume=None, problem=None, data=None):
problem = get_default(problem, self.problem)
opts = self.app_options
iw_dir = nm.array(opts.incident_wave_dir, dtype=nm.float64)
dim = problem.get_dim()
assert_(dim == iw_dir.shape[0])
iw_dir = iw_dir / nla.norm(iw_dir)
dispersion = data[self.requires[0]]
wave_vectors = dispersion.logs.eig_vectors
pas = []
iw_dir = iw_dir / nla.norm(iw_dir)
idims = list(range(iw_dir.shape[0]))
pi2 = 0.5 * nm.pi
for vecs in wave_vectors:
pa = nm.empty(vecs.shape[:-1], dtype=nm.float64)
for ir, vec in enumerate(vecs):
for ic in idims:
vv = vec[:,ic]
# Ensure the angle is in [0, pi/2].
val = nm.arccos(nm.dot(iw_dir, vv) / nla.norm(vv))
if val > pi2:
val = nm.pi - val
pa[ir,ic] = val
pas.append(pa)
return pas
[docs]class PhaseVelocity(MiniAppBase):
"""
Compute phase velocity.
"""
[docs] def process_options(self):
get = self.options.get
return Struct(eigensolver=get('eigensolver', 'eig.sgscipy'))
def __call__(self, volume=None, problem=None, data=None):
problem = get_default(problem, self.problem)
opts = self.app_options
dv_info, cat = [data[ii] for ii in self.requires]
output('average density:', dv_info.average_density)
dim = problem.get_dim()
eye = nm.eye(dim, dim, dtype=nm.float64)
mtx_mass = eye * dv_info.average_density
meigs, mvecs = eig(mtx_mass, mtx_b=cat,
eigenvectors=True, solver_kind=opts.eigensolver)
phase_velocity = 1.0 / nm.sqrt(meigs)
return phase_velocity