import numpy as nm
from sfepy.base.base import assert_
from sfepy.linalg import dot_sequences
from sfepy.terms.terms import Term, terms
from sfepy.terms.terms_dot import ScalarDotMGradScalarTerm
[docs]
class DiffusionTerm(Term):
r"""
General diffusion term with permeability :math:`K_{ij}`. Can be
evaluated. Can use derivatives.
:Definition:
.. math::
\int_{\Omega} K_{ij} \nabla_i q \nabla_j p
:Arguments:
- material: :math:`K_{ij}`
- virtual/parameter_1: :math:`q`
- state/parameter_2: :math:`p`
"""
name = 'dw_diffusion'
arg_types = (('material', 'virtual', 'state'),
('material', 'parameter_1', 'parameter_2'))
arg_shapes = {'material' : 'D, D', 'virtual' : (1, 'state'),
'state' : 1, 'parameter_1' : 1, 'parameter_2' : 1}
modes = ('weak', 'eval')
symbolic = {'expression': 'div( K * grad( u ) )',
'map' : {'u' : 'state', 'K' : 'material'}}
[docs]
def get_fargs(self, mat, virtual, state,
mode=None, term_mode=None, diff_var=None, **kwargs):
vg, _ = self.get_mapping(state)
if mat is None:
if self.name in ('dw_laplace', 'dw_st_pspg_p'):
n_el, n_qp, _, _, _ = self.get_data_shape(state)
mat = nm.ones((1, n_qp, 1, 1), dtype=nm.float64)
if mode == 'weak':
if diff_var is None:
grad = self.get(state, 'grad')
fmode = 0
else:
grad = nm.array([0], ndmin=4, dtype=nm.float64)
fmode = 1
return grad, mat, vg, fmode
elif mode == 'eval':
grad1 = self.get(virtual, 'grad')
grad2 = self.get(state, 'grad')
return grad1, grad2, mat, vg
else:
raise ValueError('unsupported evaluation mode in %s! (%s)'
% (self.name, mode))
[docs]
def get_eval_shape(self, mat, virtual, state,
mode=None, term_mode=None, diff_var=None, **kwargs):
n_el, n_qp, dim, n_en, n_c = self.get_data_shape(state)
return (n_el, 1, 1, 1), state.dtype
[docs]
def set_arg_types(self):
if self.mode == 'weak':
self.function = terms.dw_diffusion
else:
self.function = terms.d_diffusion
[docs]
class SDDiffusionTerm(Term):
r"""
Diffusion sensitivity analysis term.
:Definition:
.. math::
\int_{\Omega} \hat{K}_{ij} \nabla_i q\, \nabla_j p
.. math::
\hat{K}_{ij} = K_{ij}\left(
\delta_{ik}\delta_{jl} \nabla \cdot \ul{\Vcal}
- \delta_{ik}{\partial \Vcal_j \over \partial x_l}
- \delta_{jl}{\partial \Vcal_i \over \partial x_k}\right)
:Arguments:
- material: :math:`K_{ij}`
- parameter_q: :math:`q`
- parameter_p: :math:`p`
- parameter_mv: :math:`\ul{\Vcal}`
"""
name = 'ev_sd_diffusion'
arg_types = ('material', 'parameter_q', 'parameter_p',
'parameter_mv')
arg_shapes = {'material' : 'D, D',
'parameter_q' : 1, 'parameter_p' : 1,
'parameter_mv' : 'D'}
function = staticmethod(terms.d_sd_diffusion)
[docs]
def get_fargs(self, mat, parameter_q, parameter_p, parameter_mv,
mode=None, term_mode=None, diff_var=None, **kwargs):
vg, _ = self.get_mapping(parameter_p)
grad_q = self.get(parameter_q, 'grad')
grad_p = self.get(parameter_p, 'grad')
grad_mv = self.get(parameter_mv, 'grad')
div_mv = self.get(parameter_mv, 'div')
return grad_q, grad_p, grad_mv, div_mv, mat, vg
[docs]
def get_eval_shape(self, mat, parameter_q, parameter_p, parameter_mv,
mode=None, term_mode=None, diff_var=None, **kwargs):
n_el, n_qp, dim, n_en, n_c = self.get_data_shape(parameter_q)
return (n_el, 1, 1, 1), parameter_q.dtype
[docs]
class LaplaceTerm(DiffusionTerm):
r"""
Laplace term with :math:`c` coefficient. Can be
evaluated. Can use derivatives.
:Definition:
.. math::
\int_{\Omega} c \nabla q \cdot \nabla p
:Arguments 1:
- material: :math:`c`
- virtual/parameter_1: :math:`q`
- state/parameter_2: :math:`p`
"""
name = 'dw_laplace'
arg_types = (('opt_material', 'virtual', 'state'),
('opt_material', 'parameter_1', 'parameter_2'))
arg_shapes = [{'opt_material' : '1, 1', 'virtual' : (1, 'state'),
'state' : 1, 'parameter_1' : 1, 'parameter_2' : 1},
{'opt_material' : None}]
modes = ('weak', 'eval')
symbolic = {'expression': 'c * div( grad( u ) )',
'map' : {'u' : 'state', 'c' : 'opt_material'}}
[docs]
def set_arg_types(self):
if self.mode == 'weak':
self.function = terms.dw_laplace
else:
self.function = terms.d_laplace
[docs]
class DiffusionRTerm(Term):
r"""
Diffusion-like term with material parameter :math:`K_{j}` (to
use on the right-hand side).
:Definition:
.. math::
\int_{\Omega} K_{j} \nabla_j q
:Arguments:
- material : :math:`K_j`
- virtual : :math:`q`
"""
name = 'dw_diffusion_r'
arg_types = ('material', 'virtual')
arg_shapes = {'material' : 'D, 1', 'virtual' : (1, None)}
function = staticmethod(terms.dw_diffusion_r)
[docs]
def get_fargs(self, mat, virtual,
mode=None, term_mode=None, diff_var=None, **kwargs):
vg, _ = self.get_mapping(virtual)
return mat, vg
[docs]
class DiffusionCoupling(Term):
r"""
Diffusion copupling term with material parameter :math:`K_{j}`.
:Definition:
.. math::
\int_{\Omega} p K_{j} \nabla_j q \mbox{ , }
\int_{\Omega} q K_{j} \nabla_j p
:Arguments:
- material : :math:`K_{j}`
- virtual : :math:`q`
- state : :math:`p`
"""
name = 'dw_diffusion_coupling'
arg_types = (('material', 'virtual', 'state'),
('material', 'state', 'virtual'),
('material', 'parameter_1', 'parameter_2'))
arg_shapes = {'material' : 'D, 1', 'virtual' : (1, 'state'),
'state' : 1, 'parameter_1' : 1, 'parameter_2' : 1}
modes = ('weak0', 'weak1', 'eval')
[docs]
@staticmethod
def d_fun(out, mat, val, grad, vg):
out_qp = val * dot_sequences(mat, grad, 'ATB')
status = vg.integrate(out, out_qp)
return status
[docs]
@staticmethod
def dw_fun(out, val, mat, bf, vg, fmode):
if fmode == 0:
status = terms.mulAB_integrate(out, vg.bfg, mat * val, vg.cmap,
mode='ATB')
elif fmode == 1:
status = terms.mulAB_integrate(out, bf * mat, val, vg.cmap,
mode='ATB')
elif fmode == 2:
status = terms.mulAB_integrate(out, vg.bfg, mat * bf, vg.cmap,
mode='ATB')
elif fmode == 3:
status = terms.mulAB_integrate(out, mat * bf, vg.bfg, vg.cmap,
mode='ATB')
return status
[docs]
def get_fargs( self, mat, virtual, state,
mode=None, term_mode=None, diff_var=None, **kwargs):
vg, _ = self.get_mapping(virtual)
if mode == 'weak':
vgs, _ = self.get_mapping(state)
if diff_var is None:
if self.mode == 'weak0':
val = self.get(state, 'val')
fmode = 0
else:
val = self.get(virtual, 'grad')
fmode = 1
else:
val = nm.array([0], ndmin=4, dtype=nm.float64)
if self.mode == 'weak0':
fmode = 2
else:
fmode = 3
return val, mat, vgs.bf, vg, fmode
elif mode == 'eval':
grad = self.get(virtual, 'grad')
val = self.get(state, 'val')
return mat, val, grad, vg
else:
raise ValueError('unsupported evaluation mode in %s! (%s)'
% (self.name, mode))
[docs]
def get_eval_shape(self, mat, virtual, state,
mode=None, term_mode=None, diff_var=None, **kwargs):
n_el, n_qp, dim, n_en, n_c = self.get_data_shape(state)
return (n_el, 1, 1, 1), state.dtype
[docs]
def set_arg_types( self ):
if self.mode[:-1] == 'weak':
self.function = self.dw_fun
else:
self.function = self.d_fun
[docs]
class DiffusionVelocityTerm( Term ):
r"""
Evaluate diffusion velocity.
Supports 'eval', 'el_avg' and 'qp' evaluation modes.
:Definition:
.. math::
- \int_{\cal{D}} K_{ij} \nabla_j p
:Arguments:
- material : :math:`K_{ij}`
- parameter : :math:`p`
"""
name = 'ev_diffusion_velocity'
arg_types = ('material', 'parameter')
arg_shapes = {'material' : 'D, D', 'parameter' : 1}
integration = ('cell', 'facet_extra')
[docs]
@staticmethod
def function(out, grad, mat, vg, fmode):
dvel = dot_sequences(mat, grad)
if fmode == 2:
out[:] = dvel
status = 0
else:
status = vg.integrate(out, dvel, fmode)
out *= -1.0
return status
[docs]
def get_fargs(self, mat, parameter,
mode=None, term_mode=None, diff_var=None, **kwargs):
vg, _ = self.get_mapping(parameter)
grad = self.get(parameter, 'grad')
fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1)
return grad, mat, vg, fmode
[docs]
def get_eval_shape(self, mat, parameter,
mode=None, term_mode=None, diff_var=None, **kwargs):
n_el, n_qp, dim, n_en, n_c = self.get_data_shape(parameter)
if mode != 'qp':
n_qp = 1
return (n_el, n_qp, dim, 1), parameter.dtype
[docs]
class SurfaceFluxTerm(Term):
r"""
Surface flux term.
Supports 'eval', 'el_eval' and 'el_avg' evaluation modes.
:Definition:
.. math::
\int_{\Gamma} \ul{n} \cdot K_{ij} \nabla_j p
:Arguments:
- material: :math:`\ul{K}`
- parameter: :math:`p`,
"""
name = 'ev_surface_flux'
arg_types = ('material', 'parameter')
arg_shapes = {'material' : 'D, D', 'parameter' : 1}
integration = 'facet_extra'
function = staticmethod(terms.d_surface_flux)
[docs]
def get_fargs(self, mat, parameter,
mode=None, term_mode=None, diff_var=None, **kwargs):
sg, _ = self.get_mapping(parameter)
grad = self.get(parameter, 'grad')
fmode = {'eval' : 0, 'el_avg' : 1}.get(mode, 0)
return grad, mat, sg, fmode
[docs]
def get_eval_shape(self, mat, parameter,
mode=None, term_mode=None, diff_var=None, **kwargs):
n_fa, n_qp, dim, n_en, n_c = self.get_data_shape(parameter)
return (n_fa, 1, 1, 1), parameter.dtype
[docs]
class SurfaceFluxOperatorTerm(Term):
r"""
Surface flux operator term.
:Definition:
.. math::
\int_{\Gamma} q \ul{n} \cdot \ull{K} \cdot \nabla p
:Arguments:
- material : :math:`\ull{K}`
- virtual : :math:`q`
- state : :math:`p`
"""
name = 'dw_surface_flux'
arg_types = ('opt_material', 'virtual', 'state')
arg_shapes = [{'opt_material' : 'D, D', 'virtual' : (1, 'state'),
'state' : 1},
{'opt_material' : None}]
integration = 'facet_extra'
function = staticmethod(terms.dw_surface_flux)
[docs]
def get_fargs(self, mat, virtual, state,
mode=None, term_mode=None, diff_var=None, **kwargs):
sg, _ = self.get_mapping(state)
sd = state.field.extra_data[f'sd_{self.region.name}']
self.get_mapping(virtual) # Creates BQP for the basis.
bf = virtual.field.eval_basis(sd.bkey, 0, self.integral)
if mat is None:
_, n_qp, dim, _, _ = self.get_data_shape(state)
mat = nm.empty((1, n_qp, dim, dim), dtype=nm.float64)
mat[..., :, :] = nm.eye(dim, dtype=nm.float64)
if diff_var is None:
grad = self.get(state, 'grad', integration='facet_extra')
fmode = 0
else:
grad = nm.array([0], ndmin=4, dtype=nm.float64)
fmode = 1
return grad, mat, bf, sg, sd.fis, fmode
[docs]
class ConvectVGradSTerm(Term):
r"""
Scalar gradient term with convective velocity.
:Definition:
.. math::
\int_{\Omega} q (\ul{u} \cdot \nabla p)
:Arguments:
- virtual : :math:`q`
- state_v : :math:`\ul{u}`
- state_s : :math:`p`
"""
name = 'dw_convect_v_grad_s'
arg_types = ('virtual', 'state_v', 'state_s')
arg_shapes = [{'virtual' : (1, 'state_s'), 'state_v' : 'D', 'state_s' : 1}]
function = staticmethod(terms.dw_convect_v_grad_s)
[docs]
def get_fargs(self, virtual, state_v, state_s,
mode=None, term_mode=None, diff_var=None, **kwargs):
vvg, _ = self.get_mapping(state_v)
svg, _ = self.get_mapping(state_s)
if diff_var is None:
grad_s = self.get(state_s, 'grad')
val_v = self.get(state_v, 'val')
fmode = 0
elif diff_var == state_s.name:
grad_s = nm.array([0], ndmin=4, dtype=nm.float64)
val_v = self.get(state_v, 'val')
fmode = 1
else:
assert_(diff_var == state_v.name)
grad_s = self.get(state_s, 'grad')
val_v = nm.array([0], ndmin=4, dtype=nm.float64)
fmode = 2
return val_v, grad_s, vvg, svg, fmode
[docs]
class AdvectDivFreeTerm(ScalarDotMGradScalarTerm):
r"""
Advection of a scalar quantity :math:`p` with the advection velocity
:math:`\ul{y}` given as a material parameter (a known function of space and
time).
The advection velocity has to be divergence-free!
:Definition:
.. math::
\int_{\Omega} \nabla \cdot (\ul{y} p) q
= \int_{\Omega} (\underbrace{(\nabla \cdot \ul{y})}_{\equiv 0}
+ \ul{y} \cdot \nabla) p) q
:Arguments:
- material : :math:`\ul{y}`
- virtual : :math:`q`
- state : :math:`p`
"""
name = 'dw_advect_div_free'
arg_types = ('material', 'virtual', 'state')
arg_shapes = {'material' : 'D, 1', 'virtual' : ('1', 'state'),
'state' : '1'}
mode = 'grad_state'
[docs]
class NonlinearDiffusionTerm(Term):
r"""
The diffusion term with a scalar coefficient given by a user
supplied function of the state variable.
:Definition:
.. math::
\int_{\Omega} \nabla q \cdot \nabla p f(p)
:Arguments:
- fun : :math:`f(p)`
- dfun : :math:`\partial f(p) / \partial p`
- virtual : :math:`q`
- state : :math:`p`
"""
name = 'dw_nl_diffusion'
arg_types = ('fun', 'dfun', 'virtual', 'state')
arg_shapes = {'fun' : lambda x: x,
'dfun' : lambda x: x,
'virtual' : (1, 'state'),
'state' : 1}
[docs]
@staticmethod
def function(out, out_qp, geo):
status = geo.integrate(out, out_qp)
return status
[docs]
def get_fargs(self, fun, dfun, var1, var2,
mode=None, term_mode=None, diff_var=None, **kwargs):
vg1, _ = self.get_mapping(var1)
vg2, _ = self.get_mapping(var2)
if diff_var is None:
geo = vg1
val_grad_qp = self.get(var2, 'grad')
val_qp = fun(self.get(var2, 'val'))
out_qp = dot_sequences(vg1.bfg, val_grad_qp*val_qp,'ATB')
else:
geo = vg1
val_grad_qp = self.get(var2, 'grad')
val_d_qp = dfun(self.get(var2, 'val'))
val_qp = fun(self.get(var2, 'val'))
out_qp = (dot_sequences(vg1.bfg, vg2.bfg*val_qp,'ATB') +
dot_sequences(vg1.bfg, val_grad_qp*val_d_qp,'ATB')*vg2.bf)
return out_qp, geo