import numpy as nm
from sfepy.base.base import Struct
from sfepy.terms.terms import terms
from sfepy.terms.terms_hyperelastic_base import HyperElasticBase
_msg_missing_data = 'missing family data!'
[docs]class HyperElasticULBase(HyperElasticBase):
"""
Base class for all hyperelastic terms in UL formulation family.
The subclasses should have the following static method attributes:
- `stress_function()` (the stress)
- `tan_mod_function()` (the tangent modulus)
The common (family) data are cached in the evaluate cache of state
variable.
"""
family_function = staticmethod(terms.dq_finite_strain_ul)
weak_function = staticmethod(terms.dw_he_rtm)
fd_cache_name = 'ul_common'
hyperelastic_mode = 1
[docs] def compute_family_data(self, state):
ap, vg = self.get_approximation(state, get_saved=True)
vec = self.get_vector(state)
n_el, n_qp, dim, n_en, n_c = self.get_data_shape(state)
sym = dim * (dim + 1) / 2
shapes = {
'mtx_f' : (n_el, n_qp, dim, dim),
'det_f' : (n_el, n_qp, 1, 1),
'sym_b' : (n_el, n_qp, sym, 1),
'tr_b' : (n_el, n_qp, 1, 1),
'in2_b' : (n_el, n_qp, 1, 1),
'green_strain' : (n_el, n_qp, sym, 1),
}
data = Struct(name='ul_family_data')
for key, shape in shapes.iteritems():
setattr(data, key, nm.zeros(shape, dtype=nm.float64))
self.family_function(data.mtx_f,
data.det_f,
data.sym_b,
data.tr_b,
data.in2_b,
data.green_strain,
vec, vg, ap.econn)
return data
[docs]class NeoHookeanULTerm(HyperElasticULBase):
r"""
Hyperelastic neo-Hookean term. Effective stress :math:`\tau_{ij} = \mu
J^{-\frac{2}{3}}(b_{ij} - \frac{1}{3}b_{kk}\delta_{ij})`.
:Definition:
.. math::
\int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
:Arguments:
- material : :math:`\mu`
- virtual : :math:`\ul{v}`
- state : :math:`\ul{u}`
"""
name = 'dw_ul_he_neohook'
family_data_names = ['det_f', 'tr_b', 'sym_b']
stress_function = staticmethod(terms.dq_ul_he_stress_neohook)
tan_mod_function = staticmethod(terms.dq_ul_he_tan_mod_neohook)
# class NeoHookeanULHTerm(NeoHookeanULTerm):
# r"""
# Hyperelastic neo-Hookean term. Geometrical configuration given by
# parameter :math:`\ul{w}`. Effective stress :math:`\tau_{ij} = \mu
# J^{-\frac{2}{3}}(b_{ij} - \frac{1}{3}b_{kk}\delta_{ij})`.
# :Definition:
# .. math::
# \int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
# :Arguments 1:
# - material : :math:`\mu`
# - virtual : :math:`\ul{v}`
# - state : :math:`\ul{u}`
# - state_u : :math:`\ul{w}`
# """
# name = 'dw_ul_he_neohook_h'
# arg_types = ('material', 'virtual', 'state', 'state_u')
# use_caches = {'finite_strain_ul' : [['state_u']]}
# def __init__(self, *args, **kwargs):
# HyperElasticULBase.__init__(self, *args, **kwargs)
# self.call_mode = 1
# class NeoHookeanULEHTerm(NeoHookeanULTerm):
# r"""
# Hyperelastic neo-Hookean term.
# Geometrical configuration given by parameter :math:`\ul{w}`.
# Effective stress :math:`\tau_{ij} = \mu J^{-\frac{2}{3}}(b_{ij} - \frac{1}{3}b_{kk}\delta_{ij})`.
# :Definition:
# .. math::
# \int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
# :Arguments:
# - material : :math:`\mu`
# - parameter_1 : :math:`\ul{v}`
# - parameter_2 : :math:`\ul{u}`
# - state_u : :math:`\ul{w}`
# """
# name = 'd_ul_he_neohook_h'
# arg_types = ('material', 'parameter_1', 'parameter_2', 'state_u')
# use_caches = {'finite_strain_ul' : [['state_u']]}
# def __init__(self, *args, **kwargs):
# HyperElasticULBase.__init__(self, *args, **kwargs)
# self.call_mode = 2
# class NeoHookeanULEvalTerm(Term):
# name = 'de_ul_he_neohook'
# arg_types = ('material', 'state_u')
# use_caches = {'finite_strain_ul' : [['state_u']]}
# function = {'stress': terms.dq_ul_he_stress_neohook,
# 'element_contribution' : terms.de_he_rtm}
# def __call__( self, diff_var = None, chunk_size = None, **kwargs ):
# mat, state_u = self.get_args( ['material', 'state_u'], **kwargs )
# ap, vg = self.get_approximation(state_u)
# dim = ap.dim[0]
# sym = (dim + 1) * dim / 2
# shape = (chunk_size, 1, sym, 1)
# cache = self.get_cache('finite_strain_ul', 0)
# detF, trB, B = cache(['detF', 'trB', 'B'], self, 0, state=state_u)
# stress = nm.empty_like(B)
# fun = self.function['stress']
# fun(stress, mat, detF, trB, B)
# fun = self.function['element_contribution']
# for out, chunk in self.char_fun(chunk_size, shape):
# status = fun(out, stress, detF, vg, chunk, 1)
# out1 = nm.sum(out,0).reshape((sym,))
# yield out1, chunk, status
# class BulkPenaltyULHTerm(BulkPenaltyULTerm):
# r"""
# Hyperelastic bulk penalty term.
# Geometrical configuration given by parameter :math:`\ul{w}`.
# Stress :math:`\tau_{ij} = K(J-1)\; J \delta_{ij}`.
# :Definition:
# .. math::
# \int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
# :Arguments:
# - material : :math:`K`
# - virtual : :math:`\ul{v}`
# - state : :math:`\ul{u}`
# - state_u : :math:`\ul{w}`
# """
# name = 'dw_ul_bulk_penalty_h'
# arg_types = ('material', 'virtual', 'state', 'state_u')
# use_caches = {'finite_strain_ul' : [['state_u']]}
# def __init__(self, *args, **kwargs):
# HyperElasticULBase.__init__(self, *args, **kwargs)
# self.call_mode = 1
# class BulkPenaltyULEHTerm(BulkPenaltyULTerm):
# r"""
# Hyperelastic bulk penalty term.
# Geometrical configuration given by parameter :math:`\ul{w}`.
# Stress :math:`\tau_{ij} = K(J-1)\; J \delta_{ij}`.
# :Definition:
# .. math::
# \int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
# :Arguments:
# - material : :math:`K`
# - parameter_1 : :math:`\ul{v}`
# - parameter_2 : :math:`\ul{u}`
# - state_u : :math:`\ul{w}`
# """
# name = 'd_ul_bulk_penalty_h'
# arg_types = ('material', 'parameter_1', 'parameter_2', 'state_u')
# use_caches = {'finite_strain_ul' : [['state_u']]}
# def __init__(self, *args, **kwargs):
# HyperElasticULBase.__init__(self, *args, **kwargs)
# self.call_mode = 2
[docs]class MooneyRivlinULTerm(HyperElasticULBase):
r"""
Hyperelastic Mooney-Rivlin term.
:Definition:
.. math::
\int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
:Arguments:
- material : :math:`\kappa`
- virtual : :math:`\ul{v}`
- state : :math:`\ul{u}`
"""
name = 'dw_ul_he_mooney_rivlin'
family_data_names = ['det_f', 'tr_b', 'sym_b', 'in2_b']
stress_function = staticmethod(terms.dq_ul_he_stress_mooney_rivlin)
tan_mod_function = staticmethod(terms.dq_ul_he_tan_mod_mooney_rivlin)
[docs]class BulkPenaltyULTerm(HyperElasticULBase):
r"""
Hyperelastic bulk penalty term. Stress :math:`\tau_{ij} = K(J-1)\; J
\delta_{ij}`.
:Definition:
.. math::
\int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
:Arguments:
- material : :math:`K`
- virtual : :math:`\ul{v}`
- state : :math:`\ul{u}`
"""
name = 'dw_ul_bulk_penalty'
family_data_names = ['det_f']
stress_function = staticmethod(terms.dq_ul_he_stress_bulk)
tan_mod_function = staticmethod(terms.dq_ul_he_tan_mod_bulk)
[docs]class BulkPressureULTerm(HyperElasticULBase):
r"""
Hyperelastic bulk pressure term. Stress :math:`S_{ij} = -p J \delta_{ij}`.
:Definition:
.. math::
\int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
:Arguments:
- virtual : :math:`\ul{v}`
- state : :math:`\ul{u}`
- state_p : :math:`p`
"""
name = 'dw_ul_bulk_pressure'
arg_types = ('virtual', 'state', 'state_p')
arg_shapes = {'virtual' : ('D', 'state'), 'state' : 'D', 'state_p' : 1}
family_data_names = ['det_f', 'sym_b']
family_function = staticmethod(terms.dq_finite_strain_ul)
weak_function = staticmethod(terms.dw_he_rtm)
weak_dp_function = staticmethod(terms.dw_ul_volume)
stress_function = staticmethod(terms.dq_ul_stress_bulk_pressure)
tan_mod_u_function = staticmethod(terms.dq_ul_tan_mod_bulk_pressure_u)
[docs] def compute_data(self, family_data, mode, **kwargs):
det_f, sym_b = family_data.det_f, family_data.sym_b
p_qp = family_data.p_qp
if mode == 0:
out = nm.empty_like(sym_b)
fun = self.stress_function
elif mode == 1:
shape = list(sym_b.shape)
shape[-1] = shape[-2]
out = nm.empty(shape, dtype=nm.float64)
fun = self.tan_mod_u_function
else:
raise ValueError('bad mode! (%d)' % mode)
fun(out, p_qp, det_f)
return out
[docs] def get_fargs(self, virtual, state, state_p,
mode=None, term_mode=None, diff_var=None, **kwargs):
vgv, _ = self.get_mapping(state)
fd = self.get_family_data(state, 'ul_common', self.family_data_names)
fd.p_qp = self.get(state_p, 'val')
if mode == 'weak':
ig = self.char_fun.ig
if diff_var != state_p.name:
if diff_var is None:
stress = self.compute_data(fd, 0, **kwargs)
self.stress_cache[ig] = stress
tan_mod = nm.array([0], ndmin=4, dtype=nm.float64)
fmode = 0
else:
stress = self.stress_cache[ig]
if stress is None:
stress = self.compute_data(fd, 0, **kwargs)
tan_mod = self.compute_data(fd, 1, **kwargs)
fmode = 1
fargs = (self.weak_function,
stress, tan_mod, fd.mtx_f, fd.det_f, vgv, fmode, 1)
else:
vgs, _ = self.get_mapping(state_p)
fargs = (self.weak_dp_function, fd.det_f, vgs, vgv, 1, -1)
return fargs
elif mode == 'el_avg':
if term_mode == 'strain':
out_qp = fd.green_strain
elif term_mode == 'stress':
out_qp = self.compute_data(fd, 0, **kwargs)
else:
raise ValueError('unsupported term mode in %s! (%s)'
% (self.name, term_mode))
return self.integrate, out_qp, vgv, 1
else:
raise ValueError('unsupported evaluation mode in %s! (%s)'
% (self.name, mode))
[docs] def get_eval_shape(self, virtual, state, state_p,
mode=None, term_mode=None, diff_var=None, **kwargs):
n_el, n_qp, dim, n_en, n_c = self.get_data_shape(state)
sym = dim * (dim + 1) / 2
return (n_el, 1, sym, 1), state.dtype
# class BulkPressureULHTerm(BulkPressureULTerm):
# r"""
# Hyperelastic bulk pressure term. Stress
# :math:`S_{ij} = -p J \delta_{ij}`.
# :Definition:
# .. math::
# \int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
# :Arguments:
# - virtual : :math:`\ul{v}`
# - state : :math:`\ul{u}`
# - state_p : :math:`p`
# - state_u : :math:`w`
# """
# name = 'dw_ul_bulk_pressure_h'
# arg_types = ('virtual', 'state', 'state_p', 'state_u')
# use_caches = {'finite_strain_ul' : [['state_u']],
# 'state_in_volume_qp' : [['state_p']]}
# def __call__(self, diff_var=None, chunk_size=None, **kwargs):
# term_mode, = self.get_kwargs(['term_mode'], **kwargs)
# virtual, state, state_p, state_u = self.get_args(**kwargs)
# apv, vgv = self.get_approximation(virtual)
# aps, vgs = self.get_approximation(state_p)
# self.set_data_shape(apv, aps)
# shape, mode = self.get_shape_grad(diff_var, chunk_size)
# cache = self.get_cache('finite_strain_ul', 0)
# family_data = cache(['detF', 'B'], self, 0, state=state_u)
# ig = self.char_fun.ig
# if term_mode is None:
# if mode < 2:
# stress = self.crt_data.stress[ig]
# if stress is None:
# stress = self.compute_crt_data(family_data, 0, **kwargs)
# self.crt_data.stress[ig] = stress
# tan_mod = self.crt_data.tan_mod[ig]
# if tan_mod is None:
# tan_mod = self.compute_crt_data(family_data, 1, **kwargs)
# self.crt_data.tan_mod[ig] = tan_mod
# fun = self.function['element_contribution']
# mtxF, detF = cache(['F', 'detF'], self, 0, state=state_u)
# if mode == 0:
# vec = self.get_vector(state)
# for out, chunk in self.char_fun(chunk_size, shape):
# out2 = nm.zeros(out.shape[:-1] + (out.shape[-2],),
# dtype=nm.float64)
# status1 = fun(out2, stress, tan_mod,
# mtxF, detF, vgv, chunk, 1, 1)
# status2 = terms.he_residuum_from_mtx(out, out2, vec, apv.econn, chunk)
# yield out, chunk, status1 or status2
# else:
# for out, chunk in self.char_fun(chunk_size, shape):
# status = fun(out, stress, tan_mod,
# mtxF, detF, vgv, chunk, 1, 1)
# yield out, chunk, status
# else:
# from sfepy.base.base import debug
# debug()
# # fun = self.function['element_contribution_dp']
# # mtxF, B, detF = cache(['F', 'B', 'detF'],
# # self, 0, state=state_u)
# # bf = aps.get_base('v', 0, self.integral)
# # for out, chunk in self.char_fun(chunk_size, shape):
# # status = fun(out, bf, mtxF, B, detF, vgv, 1, chunk, 1)
# # yield -out, chunk, status
# elif term_mode == 'd_eval':
# raise NotImplementedError
# class BulkPressureULEHTerm(BulkPressureULTerm):
# r"""
# Hyperelastic bulk pressure term. Stress
# :math:`S_{ij} = -p J \delta_{ij}`.
# :Definition:
# .. math::
# \int_{\Omega} \mathcal{L}\tau_{ij}(\ul{u}) e_{ij}(\delta\ul{v})/J
# :Arguments:
# - virtual : :math:`\ul{v}`
# - state : :math:`\ul{u}`
# - state_p : :math:`p`
# - state_u : :math:`w`
# """
# name = 'd_ul_bulk_pressure_h'
# arg_types = ('virtual', 'state', 'state_p', 'state_u')
# use_caches = {'finite_strain_ul' : [['state_u']],
# 'state_in_volume_qp' : [['state_p']]}
# def __call__(self, diff_var=None, chunk_size=None, **kwargs):
# term_mode, = self.get_kwargs(['term_mode'], **kwargs)
# par1, par2, state_p, state_u = self.get_args(**kwargs)
# apv, vgv = self.get_approximation(par1)
# aps, vgs = self.get_approximation(state_p)
# self.set_data_shape(apv, aps)
# n_el, n_qp, dim, n_ep = self.data_shape_v
# shape0 = (1, dim * n_ep, dim * n_ep)
# shape = (chunk_size, 1, 1, 1)
# cache = self.get_cache('finite_strain_ul', 0)
# family_data = cache(['detF', 'B'], self, 0, state=state_u)
# ig = self.char_fun.ig
# p1 = self.get_vector(par1)
# p2 = self.get_vector(par2)
# stress = self.crt_data.stress[ig]
# if stress is None:
# stress = self.compute_crt_data(family_data, 0, **kwargs)
# self.crt_data.stress[ig] = stress
# tan_mod = self.crt_data.tan_mod[ig]
# if tan_mod is None:
# tan_mod = self.compute_crt_data(family_data, 1, **kwargs)
# self.crt_data.tan_mod[ig] = tan_mod
# fun = self.function['element_contribution']
# mtxF, detF = cache(['F', 'detF'], self, 0, state=state_u)
# for out, chunk in self.char_fun( chunk_size, shape ):
# out2 = nm.zeros((out.shape[0],) + shape0, dtype=nm.float64)
# status1 = fun(out2, stress, tan_mod,
# mtxF, detF, vgv, chunk, 1, 1)
# status2 = terms.he_eval_from_mtx(out, out2, p1, p2, apv.econn, chunk)
# out0 = nm.sum(out)
# yield out0, chunk, status1 or status2
# class BulkPressureULEvalTerm(Term):
# name = 'de_ul_bulk_pressure'
# arg_types = ('state_u', 'state_p')
# use_caches = {'finite_strain_ul' : [['state_u']],
# 'state_in_volume_qp' : [['state_p']]}
# function = {'stress': terms.dq_ul_stress_bulk_pressure,
# 'element_contribution' : terms.de_he_rtm}
# def __call__( self, diff_var = None, chunk_size = None, **kwargs ):
# state_u, state_p = self.get_args( ['state_u', 'state_p'], **kwargs )
# ap, vg = self.get_approximation(state_u)
# dim = ap.dim[0]
# sym = (dim + 1) * dim / 2
# shape = (chunk_size, 1, sym, 1)
# cache = self.get_cache( 'finite_strain_ul', 0 )
# detF, B = cache(['detF', 'B'], self, 0, state=state_u)
# cache = self.get_cache('state_in_volume_qp', 0)
# p_qp = cache('state', self, 0, state=state_p, get_vector=self.get_vector)
# stress = nm.empty_like(B)
# fun = self.function['stress']
# fun(stress, p_qp, detF)
# fun = self.function['element_contribution']
# for out, chunk in self.char_fun(chunk_size, shape):
# status = fun(out, stress, detF, vg, chunk, 1)
# out1 = nm.sum(out,0).reshape((sym,))
# yield out1, chunk, status
[docs]class VolumeULTerm(HyperElasticULBase):
r"""
Volume term (weak form) in the updated Lagrangian formulation.
:Definition:
.. math::
\begin{array}{l}
\int_{\Omega} q J(\ul{u}) \\
\mbox{volume mode: vector for } K \from \Ical_h: \int_{T_K}
J(\ul{u}) \\
\mbox{rel\_volume mode: vector for } K \from \Ical_h:
\int_{T_K} J(\ul{u}) / \int_{T_K} 1
\end{array}
:Arguments:
- virtual : :math:`q`
- state : :math:`\ul{u}`
"""
name = 'dw_ul_volume'
arg_types = ('virtual', 'state')
arg_shapes = {'virtual' : (1, None), 'state' : 'D'}
family_data_names = ['mtx_f', 'det_f']
function = staticmethod(terms.dw_ul_volume)
[docs] def get_fargs(self, virtual, state,
mode=None, term_mode=None, diff_var=None, **kwargs):
vgs, _ = self.get_mapping(virtual)
vgv, _ = self.get_mapping(state)
fd = self.get_family_data(state, 'ul_common', self.family_data_names)
if mode == 'weak':
if diff_var is None:
fmode = 0
else:
fmode = 1
elif mode == 'eval':
if term_mode == 'volume':
fmode = 2
elif term_mode == 'rel_volume':
fmode = 3
else:
raise ValueError('unsupported term evaluation mode in %s! (%s)'
% (self.name, term_mode))
else:
raise ValueError('unsupported evaluation mode in %s! (%s)'
% (self.name, mode))
return fd.det_f, vgs, vgv, 0, fmode
[docs] def get_eval_shape(self, 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]class CompressibilityULTerm(HyperElasticULBase):
r"""
Compressibility term for the updated Lagrangian formulation
:Definition:
.. math::
\int_{\Omega} 1\over \gamma p \, q
:Arguments:
- material : :math:`\gamma`
- virtual : :math:`q`
- state : :math:`p`
- parameter_u : :math:`\ul(u)`
"""
name = 'dw_ul_compressible'
arg_types = ('material', 'virtual', 'state', 'parameter_u')
arg_shapes = {'material' : '1, 1', 'virtual' : (1, 'state'), 'state' : 1,
'parameter_u' : 'D'}
family_data_names = ['mtx_f', 'det_f']
function = staticmethod(terms.dw_volume_dot_scalar)
[docs] def get_fargs(self, bulk, virtual, state, parameter_u,
mode=None, term_mode=None, diff_var=None, **kwargs):
vgp, _ = self.get_mapping(virtual)
vgs, _ = self.get_mapping(state)
vgu, _ = self.get_mapping(parameter_u)
fd = self.get_family_data(parameter_u, 'ul_common', self.family_data_names)
coef = nm.divide(bulk, fd.det_f)
if mode == 'weak':
if diff_var is None:
val_qp = self.get(state, 'val')
fmode = 0
else:
val_qp = nm.array([0], ndmin=4, dtype=nm.float64)
fmode = 1
return coef, val_qp, vgp, vgs, fmode
else:
raise ValueError('unsupported evaluation mode in %s! (%s)'
% (self.name, mode))