phononic/band_gaps_conf.py¶
Description
Configuration classes for acoustic band gaps in a strongly heterogeneous elastic body.
"""
Configuration classes for acoustic band gaps in a strongly heterogeneous
elastic body.
"""
from __future__ import absolute_import
import numpy as nm
from sfepy.base.base import get_default, import_file, Struct
from sfepy.base.conf import ProblemConf
from sfepy.discrete.fem import MeshIO
import sfepy.discrete.fem.periodic as per
from sfepy.mechanics.matcoefs import stiffness_from_lame, TransformToPlane
from sfepy.homogenization.utils import define_box_regions, get_lattice_volume
import sfepy.homogenization.coefs_base as cb
import sfepy.homogenization.coefs_phononic as cp
per.set_accuracy(1e-8)
def get_pars(dim, lam, mu):
c = stiffness_from_lame(3, lam, mu)
if dim == 2:
tr = TransformToPlane()
try:
c = tr.tensor_plane_stress(c3=c)
except:
sym = (dim + 1) * dim // 2
c = nm.zeros((sym, sym), dtype=nm.float64)
return c
def set_coef_d(variables, ir, ic, mode, pis, corrs_rs):
mode2var = {'row' : 'u1_m', 'col' : 'u2_m'}
val = pis.states[ir, ic]['u_m'] + corrs_rs.states[ir, ic]['u_m']
variables[mode2var[mode]].set_data(val)
class BandGapsConf(Struct):
"""
Configuration class for acoustic band gaps in a strongly heterogeneous
elastic body.
"""
def __init__(self, filename, approx, region_selects, mat_pars, options,
evp_options, eigenmomenta_options, band_gaps_options,
coefs_save_name='coefs',
corrs_save_names=None,
incwd=None,
output_dir=None, **kwargs):
Struct.__init__(self, approx=approx, region_selects=region_selects,
mat_pars=mat_pars, options=options,
evp_options=evp_options,
eigenmomenta_options=eigenmomenta_options,
band_gaps_options=band_gaps_options,
**kwargs)
self.incwd = get_default(incwd, lambda x: x)
self.conf = Struct()
self.conf.filename_mesh = self.incwd(filename)
output_dir = get_default(output_dir, self.incwd('output'))
default = {'evp' : 'evp', 'corrs_rs' : 'corrs_rs'}
self.corrs_save_names = get_default(corrs_save_names,
default)
io = MeshIO.any_from_filename(self.conf.filename_mesh)
self.bbox, self.dim = io.read_bounding_box(ret_dim=True)
rpc_axes = nm.eye(self.dim, dtype=nm.float64) \
* (self.bbox[1] - self.bbox[0])
self.conf.options = options
self.conf.options.update({
'output_dir' : output_dir,
'volume' : {
'value' : get_lattice_volume(rpc_axes),
},
'coefs' : 'coefs',
'requirements' : 'requirements',
'coefs_filename' : coefs_save_name,
})
self.conf.mat_pars = mat_pars
self.conf.solvers = self.define_solvers()
self.conf.regions = self.define_regions()
self.conf.materials = self.define_materials()
self.conf.fields = self.define_fields()
self.conf.variables = self.define_variables()
(self.conf.ebcs, self.conf.epbcs,
self.conf.lcbcs, self.all_periodic) = self.define_bcs()
self.conf.functions = self.define_functions()
self.conf.integrals = self.define_integrals()
self.equations, self.expr_coefs = self.define_equations()
self.conf.coefs = self.define_coefs()
self.conf.requirements = self.define_requirements()
def __call__(self):
return ProblemConf.from_dict(self.conf.__dict__,
import_file(__file__))
def define_solvers(self):
solvers = {
'ls_d' : ('ls.auto_direct', {'use_presolve' : True}),
'ls_i' : ('ls.scipy_iterative', {
'method' : 'cg',
'i_max' : 1000,
'eps_a' : 1e-12,
}),
'newton' : ('nls.newton', {
'i_max' : 1,
'eps_a' : 1e-4,
}),
}
return solvers
def define_regions(self):
regions = {
'Y' : 'all',
'Y_m' : self.region_selects.matrix,
'Y_c' : self.region_selects.inclusion,
'Gamma_mc': ('r.Y_m *v r.Y_c', 'facet'),
}
regions.update(define_box_regions(self.dim,
self.bbox[0], self.bbox[1], 1e-5))
return regions
def define_materials(self):
materials = {
'm' : ({
'D_m' : self.mat_pars.D_m,
'density_m' : self.mat_pars.density_m,
'D_c' : self.mat_pars.D_c,
'density_c' : self.mat_pars.density_c,
}, None, None, {'special_constant' : True}),
}
return materials
def define_fields(self):
fields = {
'vector_Y_m' : ('real', self.dim, 'Y_m', self.approx),
'vector_Y_c' : ('real', self.dim, 'Y_c', self.approx),
'scalar_Y' : ('real', 1, 'Y', 1),
}
return fields
def define_variables(self):
variables = {
'u_m' : ('unknown field', 'vector_Y_m'),
'v_m' : ('test field', 'vector_Y_m', 'u_m'),
'Pi' : ('parameter field', 'vector_Y_m', '(set-to-None)'),
'u1_m' : ('parameter field', 'vector_Y_m', '(set-to-None)'),
'u2_m' : ('parameter field', 'vector_Y_m', '(set-to-None)'),
'u_c' : ('unknown field', 'vector_Y_c'),
'v_c' : ('test field', 'vector_Y_c', 'u_c'),
'aux' : ('parameter field', 'scalar_Y', '(set-to-None)'),
}
return variables
def define_bcs(self):
ebcs = {
'fixed_corners' : ('Corners', {'u_m.all' : 0.0}),
'fixed_gamma_mc' : ('Gamma_mc', {'u_c.all' : 0.0}),
}
epbcs = {}
all_periodic = []
for vn in ['u_m']:
val = {'%s.all' % vn : '%s.all' % vn}
epbcs.update({
'periodic_%s_x' % vn : (['Left', 'Right'], val,
'match_y_line'),
'periodic_%s_y' % vn : (['Top', 'Bottom'], val,
'match_x_line'),
})
all_periodic.extend(['periodic_%s_x' % vn, 'periodic_%s_y' % vn])
lcbcs = {}
return ebcs, epbcs, lcbcs, all_periodic
def define_functions(self):
functions = {
'match_x_line' : (per.match_x_line,),
'match_y_line' : (per.match_y_line,),
}
return functions
def define_integrals(self):
integrals = {
'i' : 2,
}
return integrals
def define_equations(self):
equations = {}
equations['corrs_rs'] = {
'balance_of_forces' :
"""dw_lin_elastic.i.Y_m( m.D_m, v_m, u_m )
= - dw_lin_elastic.i.Y_m( m.D_m, v_m, Pi )""",
}
equations['evp'] = {
'lhs' : """dw_lin_elastic.i.Y_c( m.D_c, v_c, u_c )""",
'rhs' : """dw_dot.i.Y_c( m.density_c, v_c, u_c )""",
}
expr_coefs = {
'D' : """dw_lin_elastic.i.Y_m( m.D_m, u1_m, u2_m )""",
'VF' : """ev_volume.i.%s(aux)""",
'ema' : """ev_integrate.i.Y_c( m.density_c, u_c )""",
}
return equations, expr_coefs
def define_coefs(self):
from copy import copy
ema_options = copy(self.eigenmomenta_options)
ema_options.update({'var_name' : 'u_c'})
dispersion_options = copy(self.band_gaps_options)
dispersion_options.update({'log_save_name' : 'dispersion.log'})
coefs = {
# Basic.
'VF' : {
'regions' : ['Y_m', 'Y_c'],
'expression' : self.expr_coefs['VF'],
'class' : cb.VolumeFractions,
},
'dv_info' : {
'requires' : ['c.VF'],
'region_to_material' : {'Y_m' : ('m', 'density_m'),
'Y_c' : ('m', 'density_c'),},
'class' : cp.DensityVolumeInfo,
},
'eigenmomenta' : {
'requires' : ['evp', 'c.dv_info'],
'expression' : self.expr_coefs['ema'],
'options' : ema_options,
'class' : cp.Eigenmomenta,
},
'M' : {
'requires' : ['evp', 'c.dv_info', 'c.eigenmomenta'],
'class' : cp.AcousticMassTensor,
},
'band_gaps' : {
'requires' : ['evp', 'c.eigenmomenta', 'c.M'],
'options' : self.band_gaps_options,
'class' : cp.BandGaps,
},
# Dispersion.
'D' : {
'requires' : ['pis', 'corrs_rs'],
'expression' : self.expr_coefs['D'],
'set_variables' : set_coef_d,
'class' : cb.CoefSymSym,
},
'Gamma' : {
'requires' : ['c.D'],
'options' : {
'mode' : 'simple',
'incident_wave_dir' : None,
},
'class' : cp.ChristoffelAcousticTensor,
},
'dispersion' : {
'requires' : ['evp', 'c.eigenmomenta', 'c.M', 'c.Gamma'],
'options' : dispersion_options,
'class' : cp.BandGaps,
},
'polarization_angles' : {
'requires' : ['c.dispersion'],
'options' : {
'incident_wave_dir' : None,
},
'class' : cp.PolarizationAngles,
},
# Phase velocity.
'phase_velocity' : {
'requires' : ['c.dv_info', 'c.Gamma'],
'options' : {
'eigensolver' : 'eig.sgscipy',
},
'class' : cp.PhaseVelocity,
},
'filenames' : {},
}
return coefs
def define_requirements(self):
requirements = {
# Basic.
'evp' : {
'ebcs' : ['fixed_gamma_mc'],
'epbcs' : None,
'equations' : self.equations['evp'],
'save_name' : self.corrs_save_names['evp'],
'options' : self.evp_options,
'class' : cp.SimpleEVP,
},
# Dispersion.
'pis' : {
'variables' : ['u_m'],
'class' : cb.ShapeDimDim,
},
'corrs_rs' : {
'requires' : ['pis'],
'ebcs' : ['fixed_corners'],
'epbcs' : self.all_periodic,
'equations' : self.equations['corrs_rs'],
'set_variables' : [('Pi', 'pis', 'u_m')],
'save_name' : self.corrs_save_names['corrs_rs'],
'is_linear' : True,
'class' : cb.CorrDimDim,
'is_linear' : True,
},
}
return requirements
class BandGapsRigidConf(BandGapsConf):
"""
Configuration class for acoustic band gaps in a strongly heterogeneous
elastic body with rigid inclusions.
"""
def define_regions(self):
regions = BandGapsConf.define_regions(self)
regions['Y_cr'] = regions['Y_c']
regions.update({
'Y_r' : 'vertices by select_yr',
'Y_c' : 'r.Y_cr -c r.Y_r',
})
return regions
def define_materials(self):
materials = BandGapsConf.define_materials(self)
materials['m'][0].update({
'D_r' : self.mat_pars.D_r,
'density_r' : self.mat_pars.density_r,
})
return materials
def define_fields(self):
fields = {
'vector_Y_cr' : ('real', self.dim, 'Y_cr', self.approx),
'scalar_Y' : ('real', 1, 'Y', 1),
}
return fields
def define_variables(self):
variables = {
'u' : ('unknown field', 'vector_Y_cr'),
'v' : ('test field', 'vector_Y_cr', 'u'),
'aux' : ('parameter field', 'scalar_Y', '(set-to-None)'),
}
return variables
def define_bcs(self):
ebcs = {
'fixed_gamma_mc' : ('Gamma_mc', {'u.all' : 0.0}),
}
lcbcs ={
'rigid' : ('Y_r',{'u.all' : None}, None, 'rigid'),
}
return ebcs, {}, lcbcs, []
def define_functions(self):
functions = BandGapsConf.define_functions(self)
functions.update({
'select_yr' : (self.select_yr,),
})
return functions
def define_equations(self):
equations = {}
# dw_lin_elastic.i.Y_r( m.D_r, v, u ) should have no effect!
equations['evp'] = {
'lhs' : """dw_lin_elastic.i.Y_c( m.D_c, v, u )
+ dw_lin_elastic.i.Y_r( m.D_r, v, u )""",
'rhs' : """dw_dot.i.Y_c( m.density_c, v, u )
+ dw_dot.i.Y_r( m.density_r, v, u )""",
}
expr_coefs = {
'VF' : """ev_volume.i.%s(aux)""",
'ema' : """ev_integrate.i.Y_c( m.density_c, u )
+ ev_integrate.i.Y_r( m.density_r, u )""",
}
return equations, expr_coefs
def define_coefs(self):
from copy import copy
ema_options = copy(self.eigenmomenta_options)
ema_options.update({'var_name' : 'u'})
coefs = {
# Basic.
'VF' : {
'regions' : ['Y_m', 'Y_cr', 'Y_c', 'Y_r'],
'expression' : self.expr_coefs['VF'],
'class' : cb.VolumeFractions,
},
'dv_info' : {
'requires' : ['c.VF'],
'region_to_material' : {'Y_m' : ('m', 'density_m'),
'Y_c' : ('m', 'density_c'),
'Y_r' : ('m', 'density_r'),},
'class' : cp.DensityVolumeInfo,
},
'eigenmomenta' : {
'requires' : ['evp', 'c.dv_info'],
'expression' : self.expr_coefs['ema'],
'options' : ema_options,
'class' : cp.Eigenmomenta,
},
'M' : {
'requires' : ['evp', 'c.dv_info', 'c.eigenmomenta'],
'class' : cp.AcousticMassTensor,
},
'band_gaps' : {
'requires' : ['evp', 'c.eigenmomenta', 'c.M'],
'options' : self.band_gaps_options,
'class' : cp.BandGaps,
},
'filenames' : {},
}
return coefs
def define_requirements(self):
requirements = {
# Basic.
'evp' : {
'ebcs' : ['fixed_gamma_mc'],
'epbcs' : None,
'lcbcs' : ['rigid'],
'equations' : self.equations['evp'],
'save_name' : self.corrs_save_names['evp'],
'options' : self.evp_options,
'class' : cp.SimpleEVP,
},
}
return requirements
def clip(data, plot_range):
return nm.clip(data, *plot_range)
def clip_sqrt(data, plot_range):
return nm.clip(nm.sqrt(data), *plot_range)
def normalize(data, plot_range):
aux = nm.arctan(data)
return clip(aux, plot_range)