def get_fargs(self, mat, 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:
            if self.mode == 'grad_state':
                geo = vg1
                bf_t = vg1.bf.transpose((0, 1, 3, 2))
                val_qp = self.get(var2, 'grad')
                out_qp = bf_t * dot_sequences(mat, val_qp, 'ATB')

            else:
                geo = vg2
                val_qp = self.get(var1, 'val')
                out_qp = dot_sequences(vg2.bfg, mat, 'ATB') * val_qp

            fmode = 0

        else:
            if self.mode == 'grad_state':
                geo = vg1
                bf_t = vg1.bf.transpose((0, 1, 3, 2))
                out_qp = bf_t * dot_sequences(mat, vg2.bfg, 'ATB')

            else:
                geo = vg2
                out_qp = dot_sequences(vg2.bfg, mat, 'ATB') * vg1.bf

            fmode = 1

        return out_qp, geo, fmode

    def get_fargs(self, mat, kappa, virtual, state,
                  mode=None, term_mode=None, diff_var=None, **kwargs):
        from sfepy.discrete.variables import create_adof_conn

        geo, _ = self.get_mapping(state)

        n_fa, n_qp, dim, n_fn, n_c = self.get_data_shape(virtual)

        # Expand basis for all components.
        bf = geo.bf
        ebf = nm.zeros(bf.shape[:2] + (dim, n_fn * dim), dtype=nm.float64)
        for ir in range(dim):
            ebf[..., ir, ir*n_fn:(ir+1)*n_fn] = bf[..., 0, :]

        gmat = _build_wave_strain_op(kappa, ebf)

        if diff_var is None:
            econn = state.field.get_econn('volume', self.region)
            adc = create_adof_conn(nm.arange(state.n_dof, dtype=nm.int32),
                                   econn, n_c, 0)
            vals = state()[adc]
            # Same as nm.einsum('qij,cj->cqi', gmat[0], vals)[..., None]
            aux = dot_sequences(gmat, vals[:, None, :, None])
            out_qp = dot_sequences(gmat, dot_sequences(mat, aux), 'ATB')
            fmode = 0

        else:
            out_qp = dot_sequences(gmat, dot_sequences(mat, gmat), 'ATB')
            fmode = 1

        return out_qp, geo, fmode

    def get_fargs(self, mat, kappa, virtual, state,
                  mode=None, term_mode=None, diff_var=None, **kwargs):
        from sfepy.discrete.variables import create_adof_conn, expand_basis

        geo, _ = self.get_mapping(state)

        n_el, n_qp, dim, n_en, n_c = self.get_data_shape(virtual)

        ebf = expand_basis(geo.bf, dim)

        gmat = _build_wave_strain_op(kappa, ebf)

        if diff_var is None:
            econn = state.field.get_econn('volume', self.region)
            adc = create_adof_conn(nm.arange(state.n_dof, dtype=nm.int32),
                                   econn, n_c, 0)
            vals = state()[adc]
            # Same as nm.einsum('qij,cj->cqi', gmat[0], vals)[..., None]
            aux = dot_sequences(gmat, vals[:, None, :, None])
            out_qp = dot_sequences(gmat, dot_sequences(mat, aux), 'ATB')
            fmode = 0

        else:
            out_qp = dot_sequences(gmat, dot_sequences(mat, gmat), 'ATB')
            fmode = 1

        return out_qp, geo, fmode

    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

    def get_fargs(self, kappa, virtual, state,
                  mode=None, term_mode=None, diff_var=None, **kwargs):
        from sfepy.discrete.variables import create_adof_conn, expand_basis

        geo, _ = self.get_mapping(state)

        n_el, n_qp, dim, n_en, n_c = self.get_data_shape(virtual)

        ebf = expand_basis(geo.bf, dim)

        aux = nm.einsum('i,...ij->...j', kappa, ebf)[0, :, None, :]
        kebf = insert_strided_axis(aux, 0, n_el)

        if diff_var is None:
            econn = state.field.get_econn('volume', self.region)
            adc = create_adof_conn(nm.arange(state.n_dof, dtype=nm.int32),
                                   econn, n_c, 0)
            vals = state()[adc]
            aux = dot_sequences(kebf, vals[:, None, :, None])
            out_qp = dot_sequences(kebf, aux, 'ATB')
            fmode = 0

        else:
            out_qp = dot_sequences(kebf, kebf, 'ATB')
            fmode = 1

        return out_qp, geo, fmode

    def get_cauchy_from_2pk(self, stress_in):
        """
        Get the Cauchy stress given the second Piola-Kirchhoff stress.

        .. math::

            \sigma_{ij} = J^{-1} F_{ik} S_{kl} F_{jl}
        """
        stress_in = nm.asarray(stress_in, dtype=nm.float64)

        stress_in_full = stress_in[:,:,self.s2f,0]

        val_il = dot_sequences(self.def_grad, stress_in_full)
        val_ij = dot_sequences(val_il, self.def_grad, mode='ABT')

        stress_out_full = val_ij / self.jacobian

        sh = stress_out_full.shape
        stress_out_full.shape = (sh[0], sh[1], sh[2] * sh[3])

        self._assert_symmetry(stress_out_full)

        stress_out = nm.empty_like(stress_in)
        stress_out[...,0] = stress_out_full[:,:,self.f2s]
        return stress_out

    def d_dot(out, mat, val1_qp, val2_qp, geo):
        v1, v2 = (val1_qp, val2_qp) if val1_qp.shape[2] > 1 \
                 else (val2_qp, val1_qp)
        aux = dot_sequences(v1, mat, mode='ATB')
        vec = dot_sequences(aux, v2, mode='AB')
        status = geo.integrate(out, vec)

        return status

    def function(out, mat, vg, grad, fmode):
        bf_t = vg.bf.transpose((0, 1, 3, 2))

        if fmode == 0:
            out_qp = bf_t * dot_sequences(mat, grad, 'ATB')

        else:
            bfg = vg.bfg

            out_qp = bf_t * dot_sequences(mat, bfg, 'ATB')

        status = vg.integrate(out, out_qp)

        return status

    def get_fargs(self, mat, par_w, par_u, par_mv,
                  mode=None, term_mode=None, diff_var=None, **kwargs):
        vg, _ = self.get_mapping(par_u)

        grad_w = self.get(par_w, 'grad').transpose((0,1,3,2))
        grad_u = self.get(par_u, 'grad').transpose((0,1,3,2))
        nel, nqp, nr, nc = grad_u.shape
        strain_w = grad_w.reshape((nel, nqp, nr * nc, 1))
        strain_u = grad_u.reshape((nel, nqp, nr * nc, 1))

        mat_map = {1: nm.array([0]),
                   3: nm.array([0, 2, 2, 1]),
                   6: nm.array([0, 3, 4, 3, 1, 5, 4, 5, 2])}

        mmap = mat_map[mat.shape[-1]]
        mat_ns = mat[nm.ix_(nm.arange(nel), nm.arange(nqp),
                            mmap, mmap)]

        div_mv = self.get(par_mv, 'div')
        grad_mv = self.get(par_mv, 'grad')
        opd_mv = self.op_dv(grad_mv)

        aux = dot_sequences(mat_ns, opd_mv)
        mat_mv = mat_ns * div_mv - (aux + aux.transpose((0,1,3,2)))

        return 1.0, strain_w, strain_u, mat_mv, vg

    def __init__(self, anchor, normal, bounds):
        Struct.__init__(self, anchor=nm.array(anchor, dtype=nm.float64),
                        bounds=nm.asarray(bounds, dtype=nm.float64))
        self.normal = nm.asarray(normal, dtype=nm.float64)

        norm = nm.linalg.norm
        self.normal /= norm(self.normal)

        e3 = [0.0, 0.0, 1.0]
        dd = nm.dot(e3, self.normal)
        rot_angle = nm.arccos(dd)

        if nm.abs(rot_angle) < 1e-14:
            mtx = nm.eye(3, dtype=nm.float64)
            bounds2d = self.bounds[:, :2]

        else:
            rot_axis = nm.cross([0.0, 0.0, 1.0], self.normal)
            mtx = la.make_axis_rotation_matrix(rot_axis, rot_angle)

            mm = la.insert_strided_axis(mtx, 0, self.bounds.shape[0])
            rbounds = la.dot_sequences(mm, self.bounds)
            bounds2d = rbounds[:, :2]

        assert_(nm.allclose(nm.dot(mtx, self.normal), e3,
                            rtol=0.0, atol=1e-12))

        self.adotn = nm.dot(self.anchor, self.normal)

        self.rot_angle = rot_angle
        self.mtx = mtx
        self.bounds2d = bounds2d

def get_homog_mat(ts, coors, mode, term=None, problem=None, **kwargs):
    if problem.update_materials_flag == 2 and mode == 'qp':
        out = hyperelastic_data['homog_mat']
        return {k: nm.array(v) for k, v in six.iteritems(out)}
    elif problem.update_materials_flag == 0 or not mode == 'qp':
        return

    output('get_homog_mat')
    dim = problem.domain.mesh.dim

    update_var = problem.conf.options.mesh_update_variables[0]
    state_u = problem.equations.variables[update_var]
    state_u.field.clear_mappings()
    family_data = problem.family_data(state_u, term.region,
                                      term.integral, term.integration)

    mtx_f = family_data.mtx_f.reshape((coors.shape[0],)
                                      + family_data.mtx_f.shape[-2:])
    out = get_homog_coefs_nonlinear(ts, coors, mode, mtx_f,
                                    term=term, problem=problem,
                                    iteration=problem.iiter, **kwargs)

    out['E'] = 0.5 * (la.dot_sequences(mtx_f, mtx_f, 'ATB') - nm.eye(dim))

    hyperelastic_data['time'] = ts.step
    hyperelastic_data['homog_mat_shape'] = family_data.det_f.shape[:2]
    hyperelastic_data['homog_mat'] = \
        {k: nm.array(v) for k, v in six.iteritems(out)}

    return out

    def dw_dot(out, mat, val_qp, bfve, bfsc, geo, fmode):

        nel, nqp, dim, nc = mat.shape
        nen = bfve.shape[2]

        status1 = 0
        if fmode in [0, 1, 3]:
            aux = nm.zeros((nel, nqp, dim * nen, nc), dtype=nm.float64)
            status1 = terms.actBfT(aux, bfve, mat)

        if fmode == 0:
            status2 = terms.mulAB_integrate(out, aux, val_qp, geo, 'AB')

        if fmode == 1:
            status2 = terms.mulAB_integrate(out, aux, bfsc, geo, 'AB')

        if fmode == 2:
            aux = (bfsc * dot_sequences(mat, val_qp,
                                        mode='ATB')).transpose((0,1,3,2))
            status2 = geo.integrate(out, nm.ascontiguousarray(aux))

        if fmode == 3:
            status2 = terms.mulAB_integrate(out, bfsc, aux, geo, 'ATBT')

        return status1 and status2

    def function(out, val_qp, ebf, mat, sg, diff_var):
        normals = sg.normal
        n_fa = out.shape[0]

        ebf_t = nm.tile(ebf.transpose((0, 1, 3, 2)), (n_fa, 1, 1, 1))

        if diff_var is None:
            nu = dot_sequences(normals, val_qp, 'ATB')
            nt = dot_sequences(ebf_t, normals)
            entnu = mat * nt * nu
            status = sg.integrate(out, entnu, 0)

        else:
            nt = dot_sequences(ebf_t, normals)
            entn = mat * dot_sequences(nt, nt, 'ABT')
            status = sg.integrate(out, entn, 0)

        return status

    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')
        val_qp = dot_sequences(mat, grad, mode='ATB')

        fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1)

        return val_qp, vg, fmode

    def d_lin_prestress(self, out, strain, mat, vg, fmode):
        aux = dot_sequences(mat, strain, mode="ATB")
        if fmode == 2:
            out[:] = aux
            status = 0

        else:
            status = vg.integrate(out, aux, fmode)

        return status

    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)

        n_ep = bfg.shape[3]

        if fmode == 0:
            bts = dot_sequences(mtx_b, crt, 'ATB')

            status = geo.integrate(out, bts * h0)

            # Transform to global coordinate system, one node at
            # a time.
            for iep in range(n_ep):
                ir = slice(iep, None, n_ep)
                fn = out[:, 0, ir, 0]
                fn[:] = dot_sequences(mtx_t, fn, 'AB')


        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)

            # Transform to global coordinate system, one node at
            # a time.
            dot = dot_sequences
            for iepr in range(n_ep):
                ir = slice(iepr, None, n_ep)
                for iepc in range(n_ep):
                    ic = slice(iepc, None, n_ep)
                    fn = out[:, 0, ir, ic]
                    fn[:] = dot(dot(mtx_t, fn, 'AB'), mtx_t, 'ABT')

        return status

def transform_asm_matrices(out, mtx_t):
    """
    Transform matrix assembling contributions to global coordinate system, one
    node at a time.

    Parameters
    ----------
    out : array
        The array of matrices, transformed in-place.
    mtx_t : array
        The transposed transformation matrix :math:`T`, see
        :func:`create_transformation_matrix`.
    """
    n_ep = out.shape[-1] / mtx_t.shape[-1]
    for iepr in range(n_ep):
        ir = slice(iepr, None, n_ep)
        for iepc in range(n_ep):
            ic = slice(iepc, None, n_ep)
            fn = out[:, 0, ir, ic]
            fn[:] = dot_sequences(dot_sequences(mtx_t, fn, 'AB'), mtx_t, 'ABT')

    def d_fun(out, traction, val, sg):
        tdim = 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: # Traction vector
            aux = dot_sequences(val, traction, '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

    def function(out, val_qp, mat, vg, fmode):
        if fmode == 2:
            out[:] = dot_sequences(mat, val_qp)
            status = 0

        else:
            status = terms.de_cauchy_stress(out, val_qp, mat, vg, fmode)

        out *= -1.0

        return status

def get_tangent_stress_matrix(stress, bfg):
    """
    Get the tangent stress matrix of a thin incompressible 2D membrane
    in 3D space, given a stress.

    Parameters
    ----------
    stress : array
        The components `11, 22, 12` of the second Piola-Kirchhoff stress
        tensor, shape `(n_el, n_qp, 3, 1)`.
    bfg : array
        The in-plane base function gradients, shape `(n_el, n_qp, dim-1,
        n_ep)`.

    Returns
    -------
    mtx : array
        The tangent stress matrix, shape `(n_el, n_qp, dim*n_ep, dim*n_ep)`.
    """
    n_el, n_qp, dim, n_ep = bfg.shape
    dim += 1

    mtx = nm.zeros((n_el, n_qp, dim * n_ep, dim * n_ep), dtype=nm.float64)

    g1tg1 = dot_sequences(bfg[..., 0:1, :], bfg[..., 0:1, :], "ATB")
    g1tg2 = dot_sequences(bfg[..., 0:1, :], bfg[..., 1:2, :], "ATB")
    g2tg1 = dot_sequences(bfg[..., 1:2, :], bfg[..., 0:1, :], "ATB")
    g2tg2 = dot_sequences(bfg[..., 1:2, :], bfg[..., 1:2, :], "ATB")

    aux = (
        stress[..., 0:1, :] * g1tg1
        + stress[..., 2:3, :] * g1tg2
        + stress[..., 2:3, :] * g2tg1
        + stress[..., 1:2, :] * g2tg2
    )

    mtx[..., 0 * n_ep : 1 * n_ep, 0 * n_ep : 1 * n_ep] = aux
    mtx[..., 1 * n_ep : 2 * n_ep, 1 * n_ep : 2 * n_ep] = aux
    mtx[..., 2 * n_ep : 3 * n_ep, 2 * n_ep : 3 * n_ep] = aux

    return mtx