Source code for sfepy.terms.terms_surface

from __future__ import absolute_import
import numpy as nm

from sfepy.base.base import assert_
from sfepy.terms.terms import Term, terms
from sfepy.linalg import dot_sequences
from sfepy.mechanics.contact_bodies import ContactPlane, ContactSphere
from sfepy.mechanics.tensors import get_full_indices
from sfepy.discrete.common.extmods._geommech import geme_mulAVSB3py
from six.moves import range

##
# 22.08.2006, c
[docs] class LinearTractionTerm(Term): r""" Linear traction forces, where, depending on dimension of 'material' argument, :math:`\ull{\sigma} \cdot \ul{n}` is :math:`\bar{p} \ull{I} \cdot \ul{n}` for a given scalar pressure, :math:`\ul{f}` for a traction vector, and itself for a stress tensor. The material parameter can have one of the following shapes: 1 or (1, 1), (D, 1), (S, 1) in all modes, or (D, D) in the `eval` mode only. The symmetric tensor storage (S, 1) is as follows: in 3D S = 6 and the indices ordered as :math:`[11, 22, 33, 12, 13, 23]`, in 2D S = 3 and the indices ordered as :math:`[11, 22, 12]`. :Definition: .. math:: \int_{\Gamma} \ul{v} \cdot \ull{\sigma} \cdot \ul{n}, \int_{\Gamma} \ul{v} \cdot \ul{n}, :Arguments: - material : :math:`\ull{\sigma}` - virtual : :math:`\ul{v}` """ name = 'dw_surface_ltr' arg_types = (('opt_material', 'virtual'), ('opt_material', 'parameter')) arg_shapes = [{'opt_material' : 'S, 1', 'virtual' : ('D', None), 'parameter' : 'D'}, {'opt_material' : 'D, 1'}, {'opt_material' : '1, 1'}, {'opt_material' : 'D, D'}, {'opt_material' : None}] modes = ('weak', 'eval') integration = 'facet'
[docs] @staticmethod def d_fun(out, traction, val, sg): tdim, tdim2 = traction.shape[2:] dim = val.shape[2] sym = (dim + 1) * dim // 2 if tdim == 0: aux = dot_sequences(val, sg.normal, 'ATB') elif tdim == 1: # Pressure aux = dot_sequences(val, traction * sg.normal, 'ATB') elif tdim == dim and tdim2 == 1: # Traction vector aux = dot_sequences(val, traction, 'ATB') elif tdim == dim and tdim2 == dim: # Traction tensor - nonsymmetric trn = dot_sequences(traction, sg.normal, 'ATB') aux = dot_sequences(val, trn, 'ATB') elif tdim == sym: # Traction tensor trn, ret = geme_mulAVSB3py(traction, sg.normal) aux = dot_sequences(val, trn, 'ATB') status = sg.integrate(out, aux) return status
[docs] def get_fargs(self, traction, virtual, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(virtual) if traction is not None: n_el, _, _, _, _ = self.get_data_shape(virtual) traction = Term.tile_mat(traction, n_el) if traction is None: traction = nm.zeros((0,0,0,0), dtype=nm.float64) if mode == 'weak': return traction, sg elif mode == 'eval': val = self.get(virtual, 'val') return traction, val, sg else: raise ValueError('unsupported evaluation mode in %s! (%s)' % (self.name, mode))
[docs] def get_eval_shape(self, traction, virtual, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(virtual) return (n_el, 1, 1, 1), virtual.dtype
[docs] def set_arg_types( self ): if self.mode == 'weak': self.function = terms.dw_surface_ltr else: self.function = self.d_fun
[docs] class SDLinearTractionTerm(Term): r""" Sensitivity of the linear traction term. :Definition: .. math:: \int_{\Gamma} \ul{v} \cdot (\ull{\sigma}\, \ul{n}), \int_{\Gamma} \ul{v} \cdot \ul{n}, :Arguments: - material : :math:`\ull{\sigma}` - parameter : :math:`\ul{v}` """ name = 'ev_sd_surface_ltr' arg_types = ('opt_material', 'parameter', 'parameter_mv') arg_shapes = [{'opt_material': 'S, 1', 'parameter': 'D', 'parameter_mv': 'D'}, {'opt_material': '1, 1'}, {'opt_material': 'D, 1'}, {'opt_material': 'D, D'}, {'opt_material': None}] integration = 'facet'
[docs] @staticmethod def d_fun(out, traction, val, grad_mv, div_mv, sg): tdim, tdim2 = (None, None) if traction is None else traction.shape[2:] dim = val.shape[2] sym = (dim + 1) * dim // 2 n_el, n_qp = div_mv.shape[:2] val2 = sg.normal if tdim is None: trac = nm.tile(nm.eye(dim), (n_el, n_qp, 1, 1)) elif tdim == 1: trac = nm.tile(nm.eye(dim), (n_el, n_qp, 1, 1)) * traction elif tdim == dim and tdim2 == 1: # Traction vector trac = nm.tile(nm.eye(dim), (n_el, n_qp, 1, 1)) val2 = traction elif tdim == dim and tdim2 == dim: # Traction tensor - nonsymmetric trac = traction elif tdim == sym: # Traction tensor remap = nm.array(get_full_indices(dim)).flatten() trac = traction[..., remap, :].reshape((-1, n_qp, dim, dim)) sa_trac = trac * div_mv sa_trac -= nm.einsum('qpik,qpkj->qpij', trac, grad_mv, optimize='greedy') aux = dot_sequences(val, sa_trac, 'ATB') aux = dot_sequences(aux, val2, 'AB') status = sg.integrate(out, aux) return status
[docs] def get_fargs(self, traction, par_u, par_mv, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(par_u) val = self.get(par_u, 'val') grad_mv = self.get(par_mv, 'grad', integration='facet_extra') div_mv = self.get(par_mv, 'div', integration='facet_extra') return traction, val, grad_mv, div_mv, sg
[docs] def get_eval_shape(self, traction, par_u, par_mv, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(par_u) return (n_el, 1, 1, 1), par_u.dtype
[docs] def set_arg_types(self): self.function = self.d_fun
[docs] class ContactPlaneTerm(Term): r""" Small deformation elastic contact plane term with penetration penalty. The plane is given by an anchor point :math:`\ul{A}` and a normal :math:`\ul{n}`. The contact occurs in points that orthogonally project onto the plane into a polygon given by orthogonal projections of boundary points :math:`\{\ul{B}_i\}`, :math:`i = 1, \dots, N_B` on the plane. In such points, a penetration distance :math:`d(\ul{u}) = (\ul{X} + \ul{u} - \ul{A}, \ul{n})` is computed, and a force :math:`f(d(\ul{u})) \ul{n}` is applied. The force depends on the non-negative parameters :math:`k` (stiffness) and :math:`f_0` (force at zero penetration): - If :math:`f_0 = 0`: .. math:: f(d) = 0 \mbox{ for } d \leq 0 \;, \\ f(d) = k d \mbox{ for } d > 0 \;. - If :math:`f_0 > 0`: .. math:: f(d) = 0 \mbox{ for } d \leq -\frac{2 r_0}{k} \;, \\ f(d) = \frac{k^2}{4 r_0} d^2 + k d + r_0 \mbox{ for } -\frac{2 r_0}{k} < d \leq 0 \;, \\ f(d) = k d + f_0 \mbox{ for } d > 0 \;. In this case the dependence :math:`f(d)` is smooth, and a (small) force is applied even for (small) negative penetrations: :math:`-\frac{2 r_0}{k} < d \leq 0`. :Definition: .. math:: \int_{\Gamma} \ul{v} \cdot f(d(\ul{u})) \ul{n} :Arguments: - material_f : :math:`[k, f_0]` - material_n : :math:`\ul{n}` (special) - material_a : :math:`\ul{A}` (special) - material_b : :math:`\{\ul{B}_i\}`, :math:`i = 1, \dots, N_B` (special) - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` """ name = 'dw_contact_plane' arg_types = ('material_f', 'material_n', 'material_a', 'material_b', 'virtual', 'state') arg_shapes = {'material_f' : '1, 2', 'material_n' : '.: D', 'material_a' : '.: D', 'material_b' : '.: N, D', 'virtual' : ('D', 'state'), 'state' : 'D'} geometries = ['3_4', '3_8'] integration = 'facet' def __init__(self, *args, **kwargs): Term.__init__(self, *args, **kwargs) self.cp = None
[docs] @staticmethod def function(out, force, normal, geo, fmode): bf = geo.bf[0] nbf = bf * normal[None, :, None] nbf.shape = (bf.shape[0], bf.shape[2] * normal.shape[0]) if fmode == 0: out_qp = force * nbf[None, :, :, None] else: nbf2 = nbf[:, :, None] * nbf[:, None, :] out_qp = force * nbf2[None, :, :, :] status = geo.integrate(out, nm.ascontiguousarray(out_qp)) return status
[docs] @staticmethod def smooth_f(d, k, f0, a, eps, diff): ii = nm.where((d > eps) & (d <= 0.0))[0] if not diff: out = nm.where(d > 0.0, k * d + f0, 0.0) if len(ii): di = d[ii] out[ii] = a[ii] * di**2 + k[ii] * di + f0[ii] else: out = nm.where(d > 0.0, k, 0.0) if len(ii): out[ii] = 2 * a[ii] * d[ii] + k[ii] return out
@staticmethod def _get_force_pars(force_pars, shape): k = force_pars[..., 0] f0 = force_pars[..., 1] k.shape = f0.shape = shape ir = f0 >= 1e-14 eps = nm.where(ir, - 2.0 * f0 / k, 0.0) a = nm.zeros_like(eps) a[ir] = k[ir]**2 / (4.0 * f0[ir]) return k, f0, eps, a
[docs] def get_fargs(self, force_pars, normal, anchor, bounds, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(virtual) assert_((force_pars >= 0.0).all(), 'force parameters must be non-negative!') force_pars = Term.tile_mat(force_pars, sg.n_el) if self.cp is None: self.cp = ContactPlane(anchor, normal, bounds) qps = self.get_physical_qps() qp_coors = qps.values u_qp = self.get(state, 'val').reshape(qp_coors.shape) # Deformed QP coordinates. coors = u_qp + qp_coors force = nm.zeros(coors.shape[0], dtype=nm.float64) k, f0, eps, a = self._get_force_pars(force_pars, force.shape) # Active points in contact change with displacements! ii = self.cp.mask_points(qp_coors) if diff_var is None: if ii.any(): d = self.cp.get_distance(coors[ii]) # Force in the plane normal direction. force[ii] = self.smooth_f(d, k[ii], f0[ii], a[ii], eps[ii], 0) fmode = 0 else: if ii.any(): d = self.cp.get_distance(coors[ii]) # Force in the plane normal direction derivative. force[ii] = self.smooth_f(d, k[ii], f0[ii], a[ii], eps[ii], 1) fmode = 1 force.shape = qps.shape[:2] + (1, 1) return force, self.cp.normal, sg, fmode
[docs] class ContactSphereTerm(ContactPlaneTerm): r""" Small deformation elastic contact sphere term with penetration penalty. The sphere is given by a centre point :math:`\ul{C}` and a radius :math:`R`. The contact occurs in points that are closer to :math:`\ul{C}` than :math:`R`. In such points, a penetration distance :math:`d(\ul{u}) = R - ||\ul{X} + \ul{u} - \ul{C}||` is computed, and a force :math:`f(d(\ul{u})) \ul{n}(\ul{u})` is applied, where :math:`\ul{n}(\ul{u}) = (\ul{X} + \ul{u} - \ul{C}) / ||\ul{X} + \ul{u} - \ul{C}||`. The force depends on the non-negative parameters :math:`k` (stiffness) and :math:`f_0` (force at zero penetration): - If :math:`f_0 = 0`: .. math:: f(d) = 0 \mbox{ for } d \leq 0 \;, \\ f(d) = k d \mbox{ for } d > 0 \;. - If :math:`f_0 > 0`: .. math:: f(d) = 0 \mbox{ for } d \leq -\frac{2 r_0}{k} \;, \\ f(d) = \frac{k^2}{4 r_0} d^2 + k d + r_0 \mbox{ for } -\frac{2 r_0}{k} < d \leq 0 \;, \\ f(d) = k d + f_0 \mbox{ for } d > 0 \;. In this case the dependence :math:`f(d)` is smooth, and a (small) force is applied even for (small) negative penetrations: :math:`-\frac{2 r_0}{k} < d \leq 0`. :Definition: .. math:: \int_{\Gamma} \ul{v} \cdot f(d(\ul{u})) \ul{n}(\ul{u}) :Arguments: - material_f : :math:`[k, f_0]` - material_c : :math:`\ul{C}` (special) - material_r : :math:`R` (special) - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` """ name = 'dw_contact_sphere' arg_types = ('material_f', 'material_c', 'material_r', 'virtual', 'state') arg_shapes = {'material_f' : '1, 2', 'material_c' : '.: D', 'material_r' : '.: 1', 'virtual' : ('D', 'state'), 'state' : 'D'} geometries = ['3_4', '3_8'] integration = 'facet' def __init__(self, *args, **kwargs): Term.__init__(self, *args, **kwargs) self.cs = None
[docs] @staticmethod def function(out, force, normals, fd, geo, fmode): bf = geo.bf[0] nbf = bf * normals nbf.shape = (normals.shape[0], normals.shape[1], bf.shape[2] * normals.shape[2]) if fmode == 0: out_qp = force * nbf[..., None] else: nbf2 = nbf[..., None] * nbf[..., None, :] dim = normals.shape[2] n_ep = bf.shape[2] bb = bf[:, 0] vb = nm.zeros((bf.shape[0], dim, dim * n_ep)) for ii in range(dim): vb[:, ii, ii*n_ep:(ii+1)*n_ep] = bb ee = nm.eye(3)[None, ...] eebf2 = dot_sequences(vb, dot_sequences(ee, vb), 'ATB') out_qp = force * nbf2 if fd is not None: out_qp -= fd * (eebf2[None, :, :, :] - nbf2) status = geo.integrate(out, nm.ascontiguousarray(out_qp)) return status
[docs] def get_fargs(self, force_pars, centre, radius, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(virtual) assert_((force_pars >= 0.0).all(), 'force parameters must be non-negative!') force_pars = Term.tile_mat(force_pars, sg.n_el) if self.cs is None: self.cs = ContactSphere(centre, radius) qps = self.get_physical_qps() qp_coors = qps.values u_qp = self.get(state, 'val').reshape(qp_coors.shape) # Deformed QP coordinates. coors = u_qp + qp_coors force = nm.zeros(coors.shape[0], dtype=nm.float64) normals = nm.zeros((coors.shape[0], 3), dtype=nm.float64) fd = None k, f0, eps, a = self._get_force_pars(force_pars, force.shape) # Active points in contact change with displacements! ii = self.cs.mask_points(coors, nm.abs(eps.min())) if diff_var is None: if ii.any(): d, normals[ii] = self.cs.get_distance(coors[ii]) force[ii] = self.smooth_f(d, k[ii], f0[ii], a[ii], eps[ii], 0) fmode = 0 else: if ii.any(): d, normals[ii] = self.cs.get_distance(coors[ii]) # Force derivative. force[ii] = self.smooth_f(d, k[ii], f0[ii], a[ii], eps[ii], 1) # Force / (R - d). aux = self.smooth_f(d, k[ii], f0[ii], a[ii], eps[ii], 0) fd = nm.zeros_like(force) fd[ii] = aux / (self.cs.radius - d) fd.shape = qps.shape[:2] + (1, 1) fmode = 1 force.shape = qps.shape[:2] + (1, 1) normals.shape = qps.shape[:2] + (3, 1) return force, normals, fd, sg, fmode
[docs] class SufaceNormalDotTerm(Term): r""" "Scalar traction" term, (weak form). :Definition: .. math:: \int_{\Gamma} q \ul{c} \cdot \ul{n} :Arguments: - material : :math:`\ul{c}` - virtual : :math:`q` """ name = 'dw_surface_ndot' arg_types = (('material', 'virtual'), ('material', 'parameter')) arg_shapes = {'material' : 'D, 1', 'virtual' : (1, None), 'parameter' : 1} modes = ('weak', 'eval') integration = 'facet'
[docs] @staticmethod def dw_fun(out, material, bf, sg): bf_t = nm.tile(bf.transpose((0, 1, 3, 2)), (out.shape[0], 1, 1, 1)) bf_t = nm.ascontiguousarray(bf_t) aux = dot_sequences(material, sg.normal, 'ATB') status = sg.integrate(out, bf_t * aux) return status
[docs] @staticmethod def d_fun(out, material, val, sg): aux = dot_sequences(material, sg.normal, 'ATB') status = sg.integrate(out, val * aux) return status
[docs] def get_fargs(self, mat, virtual, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(virtual) if mode == 'weak': return mat, sg.bf, sg elif mode == 'eval': val = self.get(virtual, 'val') return mat, val, sg else: raise ValueError('unsupported evaluation mode in %s! (%s)' % (self.name, mode))
[docs] def get_eval_shape(self, mat, virtual, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(virtual) return (n_el, 1, 1, 1), virtual.dtype
[docs] def set_arg_types( self ): if self.mode == 'weak': self.function = self.dw_fun else: self.function = self.d_fun
[docs] class SDSufaceIntegrateTerm(Term): r""" Sensitivity of scalar traction. :Definition: .. math:: \int_{\Gamma} p \nabla \cdot \ul{\Vcal} :Arguments: - parameter : :math:`p` - parameter_mv : :math:`\ul{\Vcal}` """ name = 'ev_sd_surface_integrate' arg_types = ('parameter', 'parameter_mv') arg_shapes = {'parameter': 1, 'parameter_mv': 'D'} integration = 'facet'
[docs] @staticmethod def function(out, val_p, div_v, sg): status = sg.integrate(out, val_p * div_v) return status
[docs] def get_fargs(self, par, par_v, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(par) val_p = self.get(par, 'val') div_v = self.get(par_v, 'div', integration='facet_extra') return val_p, div_v, sg
[docs] def get_eval_shape(self, par, par_v, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(par_v) return (n_el, 1, 1, 1), par.dtype
[docs] class SurfaceJumpTerm(Term): r""" Interface jump condition. :Definition: .. math:: \int_{\Gamma} c\, q (p_1 - p_2) :Arguments: - material : :math:`c` - virtual : :math:`q` - state_1 : :math:`p_1` - state_2 : :math:`p_2` """ name = 'dw_jump' arg_types = ('opt_material', 'virtual', 'state_1', 'state_2') arg_shapes = [{'opt_material' : '1, 1', 'virtual' : (1, None), 'state_1' : 1, 'state_2' : 1}, {'opt_material' : None}] integration = 'facet'
[docs] @staticmethod def function(out, jump, mul, bf1, bf2, sg, fmode): bf_t = nm.tile(sg.bf.transpose((0, 1, 3, 2)), (out.shape[0], 1, 1, 1)) if fmode == 0: vec = bf_t * jump elif fmode == 1: vec = bf_t * bf1 else: vec = - bf_t * bf2 if mul is None: status = sg.integrate(out, vec) else: status = sg.integrate(out, mul * vec) return status
[docs] def get_fargs(self, coef, virtual, state1, state2, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(virtual) sg1, _ = self.get_mapping(state1) sg2, _ = self.get_mapping(state2) if diff_var is None: val1 = self.get(state1, 'val') val2 = self.get(state2, 'val') jump = val1 - val2 fmode = 0 else: jump = None if diff_var == state1.name: fmode = 1 else: fmode = 2 return jump, coef, sg1.bf, sg2.bf, sg, fmode