Source code for sfepy.terms.terms_membrane

import numpy as nm

from sfepy.base.base import assert_
from sfepy.linalg import dot_sequences
from sfepy.mechanics.tensors import dim2sym, transform_data
import sfepy.mechanics.membranes as membranes
from sfepy.terms.terms import Term

[docs] def eval_membrane_mooney_rivlin(a1, a2, mtx_c, c33, mode): """ Evaluate stress or tangent stiffness of the Mooney-Rivlin membrane. [1] Baoguo Wu, Xingwen Du and Huifeng Tan: A three-dimensional FE nonlinear analysis of membranes, Computers & Structures 59 (1996), no. 4, 601--605. """ a12 = 2.0 * a1[..., 0, 0] a22 = 2.0 * a2[..., 0, 0] sh = mtx_c.shape sym = dim2sym(sh[2]) c11 = mtx_c[..., 0, 0] c12 = mtx_c[..., 0, 1] c22 = mtx_c[..., 1, 1] pressure = c33 * (a12 + a22 * (c11 + c22)) if mode == 0: out = nm.empty((sh[0], sh[1], sym, 1)) # S_11, S_22, S_12. out[..., 0, 0] = -pressure * c22 * c33 + a12 + a22 * (c22 + c33) out[..., 1, 0] = -pressure * c11 * c33 + a12 + a22 * (c11 + c33) out[..., 2, 0] = +pressure * c12 * c33 - a22 * c12 else: out = nm.empty((sh[0], sh[1], sym, sym)) dp11 = a22 * c33 - pressure * c22 * c33 dp22 = a22 * c33 - pressure * c11 * c33 dp12 = 2.0 * pressure * c12 * c33 # D_11, D_22, D_33 out[..., 0, 0] = - 2.0 * ((a22 - pressure * c22) * c22 * c33**2 + c33 * c22 * dp11) out[..., 1, 1] = - 2.0 * ((a22 - pressure * c11) * c11 * c33**2 + c33 * c11 * dp22) out[..., 2, 2] = - a22 + pressure * (c33 + 2.0 * c12**2 * c33**2) \ + c12 * c33 * dp12 # D_21, D_31, D_32 out[..., 1, 0] = 2.0 * ((a22 - pressure * c33 - (a22 - pressure * c11) * c22 * c33**2) - c33 * c11 * dp11) out[..., 2, 0] = 2.0 * (-pressure * c12 * c22 * c33**2 + c12 * c33 * dp11) out[..., 2, 1] = 2.0 * (-pressure * c12 * c11 * c33**2 + c12 * c33 * dp22) out[..., 0, 1] = out[..., 1, 0] out[..., 0, 2] = out[..., 2, 0] out[..., 1, 2] = out[..., 2, 1] # D_12, D_13, D_23 ## out[..., 0, 1] = 2.0 * ((a22 - pressure * c33 ## - (a22 - pressure * c22) * c11 * c33**2) ## - c33 * c22 * dp22) ## out[..., 0, 2] = 2.0 * (a22 - pressure * c22) * c12 * c33**2 \ ## - c33 * c22 * dp12 ## out[..., 1, 2] = 2.0 * (a22 - pressure * c11) * c12 * c33**2 \ ## - c33 * c11 * dp12 return out
[docs] class TLMembraneTerm(Term): r""" Mooney-Rivlin membrane with plain stress assumption. The membrane has a uniform initial thickness :math:`h_0` and obeys a hyperelastic material law with strain energy by Mooney-Rivlin: :math:`\Psi = a_1 (I_1 - 3) + a_2 (I_2 - 3)`. :Arguments: - material_a1 : :math:`a_1` - material_a2 : :math:`a_2` - material_h0 : :math:`h_0` - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` """ name = 'dw_tl_membrane' arg_types = ('material_a1', 'material_a2', 'material_h0', 'virtual', 'state') arg_shapes = {'material_a1' : '1, 1', 'material_a2' : '1, 1', 'material_h0' : '1, 1', 'virtual' : ('D', 'state'), 'state' : 'D'} geometries = ['3_4', '3_8'] integration = 'facet'
[docs] @staticmethod def function(out, fun, *args): """ Notes ----- `fun` is either `weak_function` or `eval_function` according to evaluation mode. """ return fun(out, *args)
[docs] @staticmethod def weak_function(out, a1, a2, h0, mtx_c, c33, mtx_b, mtx_t, bfg, geo, fmode): crt = eval_membrane_mooney_rivlin(a1, a2, mtx_c, c33, fmode) if fmode == 0: bts = dot_sequences(mtx_b, crt, 'ATB') status = geo.integrate(out, bts * h0) membranes.transform_asm_vectors(out, mtx_t) else: btd = dot_sequences(mtx_b, crt, 'ATB') btdb = dot_sequences(btd, mtx_b) stress = eval_membrane_mooney_rivlin(a1, a2, mtx_c, c33, 0) kts = membranes.get_tangent_stress_matrix(stress, bfg) mtx_k = kts + btdb status = geo.integrate(out, mtx_k * h0) membranes.transform_asm_matrices(out, mtx_t) return status
[docs] @staticmethod def eval_function(out, a1, a2, h0, mtx_c, c33, mtx_b, mtx_t, geo, term_mode, fmode): if term_mode == 'strain': out_qp = membranes.get_green_strain_sym3d(mtx_c, c33) elif term_mode == 'stress': n_el, n_qp, dm, _ = mtx_c.shape dim = dm + 1 sym = dim2sym(dim) out_qp = nm.zeros((n_el, n_qp, sym, 1), dtype=mtx_c.dtype) stress = eval_membrane_mooney_rivlin(a1, a2, mtx_c, c33, 0) out_qp[..., 0:2, 0] = stress[..., 0:2, 0] out_qp[..., 3, 0] = stress[..., 2, 0] status = geo.integrate(out, out_qp, fmode) out[:, 0, :, 0] = transform_data(out.squeeze(), mtx=mtx_t) return status
def __init__(self, *args, **kwargs): Term.__init__(self, *args, **kwargs) self.mtx_t = None self.membrane_geo = None self.bfg = None
[docs] def get_fargs(self, a1, a2, h0, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): vv, vu = virtual, state sg, _ = self.get_mapping(vv) sd = vv.field.extra_data[f'sd_{self.region.name}'] if self.mtx_t is None: aux = membranes.describe_geometry(vu.field, self.region, self.integral) self.mtx_t, self.membrane_geo = aux # Transformed base function gradient w.r.t. material coordinates # in quadrature points. self.bfg = self.membrane_geo.bfg mtx_t = self.mtx_t bfg = self.bfg geo = self.membrane_geo # Displacements of element nodes. vec_u = vu.get_state_in_region(self.region) el_u = vec_u[sd.leconn] # Transform displacements to the local coordinate system. # u_new = T^T u el_u_loc = dot_sequences(el_u, mtx_t, 'AB') ## print el_u_loc mtx_c, c33, mtx_b = membranes.describe_deformation(el_u_loc, bfg) if mode == 'weak': fmode = diff_var is not None return (self.weak_function, a1, a2, h0, mtx_c, c33, mtx_b, mtx_t, bfg, geo, fmode) else: fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1) assert_(term_mode in ['strain', 'stress']) return (self.eval_function, a1, a2, h0, mtx_c, c33, mtx_b, mtx_t, geo, term_mode, fmode)
[docs] def get_eval_shape(self, a1, a2, h0, 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) sym = dim2sym(dim) return (n_el, 1, sym, 1), state.dtype