import numpy as nm

from sfepy.linalg import dot_sequences
from sfepy.terms.terms import Term, terms
from sfepy.terms.terms_th import THTerm, ETHTerm
from sfepy.terms.terms_elastic import CauchyStressTerm

class BiotTerm(Term):
    r"""
    Biot coupling term with :math:`\alpha_{ij}` given in:

    * vector form exploiting symmetry - in 3D it has the
      indices ordered as :math:`[11, 22, 33, 12, 13, 23]`, in 2D it has
      the indices ordered as :math:`[11, 22, 12]`,

    * matrix form - non-symmetric coupling parameter.

    Corresponds to weak forms of Biot gradient and divergence terms.
    Can be evaluated. Can use derivatives.

    :Definition:

    .. math::
        \int_{\Omega}  p\ \alpha_{ij} e_{ij}(\ul{v}) \mbox{ , } \int_{\Omega}
        q\ \alpha_{ij} e_{ij}(\ul{u})

    :Arguments 1:
        - material : :math:`\alpha_{ij}`
        - virtual  : :math:`\ul{v}`
        - state    : :math:`p`

    :Arguments 2:
        - material : :math:`\alpha_{ij}`
        - state    : :math:`\ul{u}`
        - virtual  : :math:`q`

    :Arguments 3:
        - material    : :math:`\alpha_{ij}`
        - parameter_v : :math:`\ul{u}`
        - parameter_s : :math:`p`
    """
    name = 'dw_biot'
    arg_types = (('material', 'virtual', 'state'),
                 ('material', 'state', 'virtual'),
                 ('material', 'parameter_v', 'parameter_s'))
    arg_shapes = [{'material' : 'S, 1',
                  'virtual/grad' : ('D', None), 'state/grad' : 1,
                  'virtual/div' : (1, None), 'state/div' : 'D',
                  'parameter_v' : 'D', 'parameter_s' : 1},
                  {'material' : 'D, D'}]

    modes = ('grad', 'div', 'eval')

    def get_fargs(self, mat, vvar, svar,
                  mode=None, term_mode=None, diff_var=None, **kwargs):

        sym_mode = False if mat.shape[-2] == mat.shape[-1] > 1 else True
        if not sym_mode:
            sh = mat.shape
            # the gradient given by 'self.get' is transposed
            mat = nm.swapaxes(mat, 2, 3)
            mat = mat.reshape(sh[:2] + (sh[2]**2, 1))

        if self.mode == 'grad':
            qp_var, qp_name = svar, 'val'

        else:
            if sym_mode:
                qp_var, qp_name = vvar, 'cauchy_strain'
            else:
                qp_var, qp_name = vvar, 'grad'

        if mode == 'weak':
            vvg, _ = self.get_mapping(vvar)
            svg, _ = self.get_mapping(svar)

            if diff_var is None:
                val_qp = self.get(qp_var, qp_name)
                if qp_name == 'grad':
                    sh = val_qp.shape
                    val_qp = val_qp.reshape(sh[:2] + (sh[2]**2, 1))

                fmode = 0

            else:
                val_qp = nm.array([0], ndmin=4, dtype=nm.float64)
                fmode = 1

            return 1.0, val_qp, mat, svg, vvg, fmode

        elif mode == 'eval':
            vvg, _ = self.get_mapping(vvar)

            if sym_mode:
                strain = self.get(vvar, 'cauchy_strain')
            else:
                strain = self.get(vvar, 'grad')
                sh = strain.shape
                strain = strain.reshape(sh[:2] + (sh[2]**2, 1))

            pval = self.get(svar, 'val')

            return 1.0, pval, strain, mat, vvg

        else:
            raise ValueError('unsupported evaluation mode in %s! (%s)'
                             % (self.name, mode))

    def get_eval_shape(self, mat, vvar, svar,
                       mode=None, term_mode=None, diff_var=None, **kwargs):
        n_el, n_qp, dim, n_en, n_c = self.get_data_shape(vvar)

        return (n_el, 1, 1, 1), vvar.dtype

    def set_arg_types(self):
        self.function = {
            'grad' : terms.dw_biot_grad,
            'div' : terms.dw_biot_div,
            'eval' : terms.d_biot_div,
        }[self.mode]

class BiotStressTerm(CauchyStressTerm):
    r"""
    Evaluate Biot stress tensor.

    It is given in the usual vector form exploiting symmetry: in 3D it has 6
    components with the indices ordered as :math:`[11, 22, 33, 12, 13, 23]`, in
    2D it has 3 components with the indices ordered as :math:`[11, 22, 12]`.

    Supports 'eval', 'el_avg' and 'qp' evaluation modes.

    :Definition:

    .. math::
        - \int_{\Omega} \alpha_{ij} \bar{p}

    .. math::
        \mbox{vector for } K \from \Ical_h:
        - \int_{T_K} \alpha_{ij} \bar{p} / \int_{T_K} 1

    .. math::
        - \alpha_{ij} \bar{p}|_{qp}

    :Arguments:
        - material  : :math:`\alpha_{ij}`
        - parameter : :math:`\bar{p}`
    """
    name = 'ev_biot_stress'
    arg_types = ('material', 'parameter')
    arg_shapes = {'material' : 'S, 1', 'parameter' : 1}

    @staticmethod
    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)