Source code for sfepy.terms.terms_new

"""
todo:
    - get row variable, col variable (if diff_var)
    - determine out shape
    - set current group to all variable arguments
    - loop over row/col dofs:
        - call term

? how to deal with components of (vector) variables?

(*) for given group, Variable has to be able to:
    - evaluate value in quadrature points
    - evaluate gradient in quadrature points
    - ? evaluate divergence in quadrature points
    - ? lazy evaluation, cache!

?? base function gradients in space elements stored now in terms - in
geometry, shared dict of geometries belongs to Equations
-> where to cache stuff? - in variables!
"""
import numpy as nm

from sfepy.base.base import output
from sfepy.terms.terms import Term, get_shape_kind
from sfepy.terms.utils import get_range_indices
from sfepy.mechanics.tensors import get_full_indices
from sfepy.linalg import dot_sequences as dot

[docs]class NewTerm(Term):
[docs] def get_geometry_key(self, variable): is_trace = self.arg_traces[variable.name] geometry_type = self.geometry_types[variable.name] region_name, iorder, ig = self.get_current_group() if is_trace: region, ig_map, ig_map_i = self.region.get_mirror_region() region_name = region.name ig = ig_map_i[ig] ap = variable.get_approximation(ig) key = (region_name, iorder, geometry_type, ap.name) return key, ig
[docs] def get_geometry(self, variable): key, ig = self.get_geometry_key(variable) geo = self.get_mapping(variable)[0] return geo, key, ig
[docs] def set_current_group(self, ig): """ Set current group for the term and all variables in its arguments. """ self.char_fun.set_current_group(ig) shape_kind = get_shape_kind(self.integration) for var in self.get_variables(): geo, geo_key, geo_ig = self.get_geometry(var) var.setup_bases(geo_key, geo_ig, geo, self.integral, shape_kind) var.set_current_group(geo_key, geo_ig)
[docs] def integrate(self, val_qp, variable): shape_kind = get_shape_kind(self.integration) geo, _, _ = self.get_geometry(variable) sh = val_qp.shape val = nm.zeros((sh[0], 1, sh[2], sh[3]), dtype=val_qp.dtype) if shape_kind == 'volume': geo.integrate(val, val_qp) else: geo.integrate(val, val_qp) return val
[docs] def evaluate(self, mode='eval', diff_var=None, **kwargs): shape_kind = get_shape_kind(self.integration) if mode == 'eval': var = self.get_variables()[0] val = 0.0 for ig in self.iter_groups(): args = self.get_args(**kwargs) val_qp = self(*args, **kwargs) _val = self.integrate(val_qp, var) val += self.sign * _val.sum() elif mode in ('el_avg', 'qp'): raise NotImplementedError() elif mode == 'weak': varr = self.get_virtual_variable() vals = [] iels = [] if diff_var is None: for ig in self.iter_groups(): args = self.get_args(**kwargs) aux = varr.get_data_shape(ig, self.integral, shape_kind, self.region.name) n_elr, n_qpr, dim, n_enr, n_cr = aux n_row = n_cr * n_enr shape = (n_elr, 1, n_row, 1) val = nm.zeros(shape, dtype=varr.dtype) for ir in varr.iter_dofs(): irs = slice(ir, ir + 1) try: val_qp = self(*args, **kwargs) except ValueError: output('%s term evaluation failed!' % self.name) raise _val = self.integrate(val_qp, varr) val[..., irs, :] = _val vals.append(self.sign * val) iels.append((ig, nm.arange(n_elr, dtype=nm.int32))) else: varc = self.get_variables(as_list=False)[diff_var] for ig in self.iter_groups(): args = self.get_args(**kwargs) aux = varr.get_data_shape(ig, self.integral, shape_kind, self.region.name) n_elr, n_qpr, dim, n_enr, n_cr = aux n_row = n_cr * n_enr aux = varc.get_data_shape(ig, self.integral, shape_kind, self.region.name) n_elc, n_qpc, dim, n_enc, n_cc = aux n_col = n_cc * n_enc shape = (n_elr, 1, n_row, n_col) val = nm.zeros(shape, dtype=varr.dtype) for ir in varr.iter_dofs(): irs = slice(ir, ir + 1) for ic in varc.iter_dofs(): ics = slice(ic, ic + 1) try: val_qp = self(*args, **kwargs) except ValueError: output('%s term evaluation failed!' % self.name) raise _val = self.integrate(val_qp, varr) val[..., irs, ics] = _val vals.append(self.sign * val) iels.append((ig, nm.arange(n_elr, dtype=nm.int32))) # Setup return value. if mode == 'eval': out = (val,) else: out = (vals, iels) # Hack: add zero status. out = out + (0,) if len(out) == 1: out = out[0] return out
[docs]class NewDiffusionTerm(NewTerm): """ """ name = 'dw_new_diffusion' arg_types = ('material', 'virtual', 'state') def __call__(self, mat, virtual, state, **kwargs): val = dot(virtual.grad(), dot(mat, state.grad()), 'ATB') return val
[docs]class NewMassScalarTerm(NewTerm): """ """ name = 'dw_new_mass_scalar' arg_types = ('virtual', 'state') def __call__(self, virtual, state, **kwargs): val = virtual.val() * state.val() return val
[docs]class NewMassTerm(NewTerm): """ Works for both scalar and vector variables. """ name = 'dw_new_mass' arg_types = ('virtual', 'state') def __call__(self, virtual, state, **kwargs): rindx = virtual.get_component_indices() cindx = state.get_component_indices() val = virtual.get_element_zeros() for ir, irs in rindx: for ic, ics in cindx: if ir == ic: val += virtual.val(ir) * state.val(ic) return val
[docs]class NewLinearElasticTerm(NewTerm): """ """ name = 'dw_new_lin_elastic' arg_types = ('material', 'virtual', 'state') def __call__(self, mat, virtual, state, **kwargs): """ Doubled out-of-diagonal strain entries! """ rindx = virtual.get_component_indices() cindx = state.get_component_indices() kindx = lindx = get_range_indices(state.dim) fi = nm.array(get_full_indices(state.dim)) val = virtual.get_element_zeros() for ir, irs in rindx: for ik, iks in kindx: irk = fi[ir, ik] irks = slice(irk, irk + 1) erk = virtual.grad(ir, ik) for ic, ics in cindx: for il, ils in lindx: icl = fi[ic, il] icls = slice(icl, icl + 1) ecl = state.grad(ic, il) val += mat[..., irks, icls] * erk * ecl return val