"""
Composite Laminate Module (:mod:`compmech.composite.laminate`)
==============================================================

.. currentmodule:: compmech.composite.laminate

"""
from __future__ import division, absolute_import

import numpy as np

from .lamina import Lamina
from .matlamina import read_laminaprop
from compmech.constants import DOUBLE
from compmech.logger import *


def read_stack(stack, plyt=None, laminaprop=None, plyts=[], laminaprops=[],
               offset=0.):
    """Read a laminate stacking sequence data.

    An ``Laminate`` object is returned based on the inputs given.

    Parameters
    ----------
    stack : list
        Angles of the stacking sequence in degrees.
    plyt : float, optional
        When all plies have the same thickness, ``plyt`` can be supplied.
    laminaprop : tuple, optional
        When all plies have the same material properties, ``laminaprop``
        can be supplied.
    plyts : list, optional
        A list of floats with the thickness of each ply.
    laminaprops : list, optional
        A list of tuples with a laminaprop for each ply.
    offset : float, optional
        Offset along the normal axis about the mid-surface, which influences
        the laminate properties.

    Notes
    -----
    ``plyt`` or ``plyts`` must be supplied
    ``laminaprop`` or ``laminaprops`` must be supplied

    For orthotropic plies, the ``laminaprop`` should be::

        laminaprop = (E11, E22, nu12, G12, G13, G23)

    For isotropic pliey, the ``laminaprop`` should be::

        laminaprop = (E, E, nu)

    """
    lam = Laminate()
    lam.offset = offset
    lam.stack = stack

    if not plyts:
        if not plyt:
            error('plyt or plyts must be supplied')
            raise ValueError
        else:
            plyts = [plyt for i in stack]

    if not laminaprops:
        if not laminaprop:
            error('laminaprop or laminaprops must be supplied')
            raise ValueError
        else:
            laminaprops = [laminaprop for i in stack]

    lam.plies = []
    for plyt, laminaprop, theta in zip(plyts, laminaprops, stack):
        laminaprop = laminaprop
        ply = Lamina()
        ply.theta = float(theta)
        ply.t = plyt
        ply.matobj = read_laminaprop(laminaprop)
        lam.plies.append(ply)

    lam.rebuild()
    lam.calc_constitutive_matrix()

    return lam


def read_lamination_parameters(thickness, laminaprop,
                               xiA1, xiA2, xiA3, xiA4,
                               xiB1, xiB2, xiB3, xiB4,
                               xiD1, xiD2, xiD3, xiD4,
                               xiE1, xiE2, xiE3, xiE4):
    r"""Calculates a laminate based on the lamination parameters.

    The lamination parameters:
    `\xi_{A1} \cdots \xi_{A4}`,  `\xi_{B1} \cdots \xi_{B4}`,
    `\xi_{C1} \cdots \xi_{C4}`,  `\xi_{D1} \cdots \xi_{D4}`,
    `\xi_{E1} \cdots \xi_{E4}`

    are used to calculate the laminate constitutive matrix.

    Parameters
    ----------
    thickness : float
        The total thickness of the laminate
    laminaprop : tuple
        The laminaprop tuple used to define the laminate material.
    xiA1 to xiD4 : float
        The 16 lamination parameters used to define the laminate.

    Returns
    -------
    lam : Laminate
        laminate with the ABD and ABDE matrices already calculated

    """
    lam = Laminate()
    lam.t = thickness
    lam.matobj = read_laminaprop(laminaprop)
    lam.xiA = np.array([1, xiA1, xiA2, xiA3, xiA4], dtype=DOUBLE)
    lam.xiB = np.array([0, xiB1, xiB2, xiB3, xiB4], dtype=DOUBLE)
    lam.xiD = np.array([1, xiD1, xiD2, xiD3, xiD4], dtype=DOUBLE)
    lam.xiE = np.array([1, xiE1, xiE2, xiE3, xiE4], dtype=DOUBLE)

    lam.calc_ABDE_from_lamination_parameters()
    return lam


class Laminate(object):
    r"""
    =========  ===========================================================
    attribute  description
    =========  ===========================================================
    plies      list of plies
    t          total thickness of the laminate
    offset     offset at the normal direction
    e1         equivalent laminate modulus in 1 direction
    e2         equivalent laminate modulus in 2 direction
    g12        equivalent laminate shear modulus in 12 direction
    nu12       equivalent laminate Poisson ratio in 12 direction
    nu21       equivalent laminate Poisson ratio in 21 direction
    xiA        laminate parameters for extensional matrix A
    xiB        laminate parameters for extension-bending matrix B
    xiD        laminate parameters for bending matrix D
    A          laminate extension matrix
    B          laminate extension-bending matrix
    D          laminate bending matrix
    E          laminate transferse shear matrix
    ABD        laminate ABD matrix
    ABDE       laminate ABD matrix with transverse shear terms
    =========  ===========================================================

    """
    def __init__(self):
        self.plies = []
        self.matobj = None
        self.t = None
        self.offset = 0.
        self.e1 = None
        self.e2 = None
        self.e3 = None
        self.nu12 = None
        self.g12 = None
        self.g13 = None
        self.g23 = None
        self.xiA = None
        self.xiB = None
        self.xiD = None
        self.A = None
        self.B = None
        self.D = None
        self.E = None
        self.ABD = None
        self.ABDE = None


    def rebuild(self):
        lam_thick = 0
        for ply in self.plies:
            ply.rebuild()
            lam_thick += ply.t
        self.t = lam_thick


    def calc_equivalent_modulus(self):
        """Calculates the equivalent laminate properties.

        The following attributes are calculated:
            e1, e2, g12, nu12, nu21

        """
        AI = np.matrix(self.ABD, dtype=DOUBLE).I
        a11, a12, a22, a33 = AI[0,0], AI[0,1], AI[1,1], AI[2,2]
        self.e1 = 1./(self.t*a11)
        self.e2 = 1./(self.t*a22)
        self.g12 = 1./(self.t*a33)
        self.nu12 = - a12 / a11
        self.nu21 = - a12 / a22


    def calc_lamination_parameters(self):
        """Calculate the lamination parameters.

        The following attributes are calculated:
            xiA, xiB, xiD, xiE

        """
        xiA1, xiA2, xiA3, xiA4 = 0, 0, 0, 0
        xiB1, xiB2, xiB3, xiB4 = 0, 0, 0, 0
        xiD1, xiD2, xiD3, xiD4 = 0, 0, 0, 0
        xiE1, xiE2, xiE3, xiE4 = 0, 0, 0, 0

        lam_thick = sum([ply.t for ply in self.plies])
        self.t = lam_thick

        h0 = -lam_thick/2. + self.offset
        for ply in self.plies:
            hk_1 = h0
            h0 += ply.t
            hk = h0

            Afac = ply.t / lam_thick
            Bfac = (2. / lam_thick**2) * (hk**2 - hk_1**2)
            Dfac = (4. / lam_thick**3) * (hk**3 - hk_1**3)
            Efac = (1. / lam_thick) * (hk - hk_1)# * (5./6) * (5./6)

            cos2t = ply.cos2t
            cos4t = ply.cos4t
            sin2t = ply.sin2t
            sin4t = ply.sin4t

            xiA1  += Afac * cos2t
            xiA2  += Afac * sin2t
            xiA3  += Afac * cos4t
            xiA4  += Afac * sin4t

            xiB1  += Bfac * cos2t
            xiB2  += Bfac * sin2t
            xiB3  += Bfac * cos4t
            xiB4  += Bfac * sin4t

            xiD1  += Dfac * cos2t
            xiD2  += Dfac * sin2t
            xiD3  += Dfac * cos4t
            xiD4  += Dfac * sin4t

            xiE1  += Efac * cos2t
            xiE2  += Efac * sin2t
            xiE3  += Efac * cos4t
            xiE4  += Efac * sin4t

        self.xiA = np.array([1, xiA1, xiA2, xiA3, xiA4], dtype=DOUBLE)
        self.xiB = np.array([0, xiB1, xiB2, xiB3, xiB4], dtype=DOUBLE)
        self.xiD = np.array([1, xiD1, xiD2, xiD3, xiD4], dtype=DOUBLE)
        self.xiE = np.array([1, xiE1, xiE2, xiE3, xiE4], dtype=DOUBLE)


    def calc_ABDE_from_lamination_parameters(self):
        """Use the ABDE matrix based on lamination parameters.

        Given the lamination parameters ``xiA``, ``xiB``, ``xiC`` and ``xiD``,
        the ABD matrix is calculated.

        """
        # dummies used to unpack vector results
        du1, du2, du3, du4, du5, du6 = 0, 0, 0, 0, 0, 0
        # A matrix terms
        A11,A22,A12, du1,du2,du3, A66,A16,A26 =\
            (self.t       ) * np.dot(self.matobj.u, self.xiA)
        # B matrix terms
        B11,B22,B12, du1,du2,du3, B66,B16,B26 =\
            (self.t**2/4. ) * np.dot(self.matobj.u, self.xiB)
        # D matrix terms
        D11,D22,D12, du1,du2,du3, D66,D16,D26 =\
            (self.t**3/12.) * np.dot(self.matobj.u, self.xiD)
        # E matrix terms
        du1,du2,du3, E44,E55,E45, du4,du5,du6 =\
            (self.t       ) * np.dot(self.matobj.u, self.xiE)

        self.A = np.array([[A11, A12, A16],
                           [A12, A22, A26],
                           [A16, A26, A66]], dtype=DOUBLE)

        self.B = np.array([[B11, B12, B16],
                           [B12, B22, B26],
                           [B16, B26, B66]], dtype=DOUBLE)

        self.D = np.array([[D11, D12, D16],
                           [D12, D22, D26],
                           [D16, D26, D66]], dtype=DOUBLE)

        # printing E acoordingly to Reddy definition for E44, E45 and E55
        self.E = np.array([[E55, E45],
                           [E45, E44]], dtype=DOUBLE)

        self.ABD = np.array([[A11, A12, A16, B11, B12, B16],
                             [A12, A22, A26, B12, B22, B26],
                             [A16, A26, A66, B16, B26, B66],
                             [B11, B12, B16, D11, D12, D16],
                             [B12, B22, B26, D12, D22, D26],
                             [B16, B26, B66, D16, D26, D66]], dtype=DOUBLE)

        # printing ABDE acoordingly to Reddy definition for E44, E45 and E55
        self.ABDE = np.array([[A11, A12, A16, B11, B12, B16, 0, 0],
                              [A12, A22, A26, B12, B22, B26, 0, 0],
                              [A16, A26, A66, B16, B26, B66, 0, 0],
                              [B11, B12, B16, D11, D12, D16, 0, 0],
                              [B12, B22, B26, D12, D22, D26, 0, 0],
                              [B16, B26, B66, D16, D26, D66, 0, 0],
                              [0, 0, 0, 0, 0, 0, E55, E45],
                              [0, 0, 0, 0, 0, 0, E45, E44]],
                               dtype=DOUBLE)


    def calc_constitutive_matrix(self):
        """Calculates the laminate constitutive matrix

        This is the commonly called ``ABD`` matrix with ``shape=(6, 6)`` when
        the classical laminated plate theory is used, or the ``ABDE`` matrix
        when the first-order shear deformation theory is used, containing the
        transverse shear terms.

        """
        self.A_general = np.zeros([5,5], dtype=DOUBLE)
        self.B_general = np.zeros([5,5], dtype=DOUBLE)
        self.D_general = np.zeros([5,5], dtype=DOUBLE)

        lam_thick = sum([ply.t for ply in self.plies])
        self.t = lam_thick

        h0 = -lam_thick/2 + self.offset
        for ply in self.plies:
            hk_1 = h0
            h0 += ply.t
            hk = h0
            self.A_general += ply.QL*(hk - hk_1)
            self.B_general += 1/2.*ply.QL*(hk**2 - hk_1**2)
            self.D_general += 1/3.*ply.QL*(hk**3 - hk_1**3)

        self.A = self.A_general[0:3, 0:3]
        self.B = self.B_general[0:3, 0:3]
        self.D = self.D_general[0:3, 0:3]
        self.E = self.A_general[3:5, 3:5]

        conc1 = np.concatenate([self.A, self.B], axis=1)
        conc2 = np.concatenate([self.B, self.D], axis=1)

        self.ABD = np.concatenate([conc1, conc2], axis=0)
        self.ABDE = np.zeros((8, 8), dtype=DOUBLE)
        self.ABDE[0:6, 0:6] = self.ABD
        self.ABDE[6:8, 6:8] = self.E


    def force_balanced_LP(self):
        r"""Force balanced lamination parameters

        The lamination parameters `\xi_{A2}` and `\xi_{A4}` are set to null
        to force a balanced laminate.

        """
        dummy, xiA1, xiA2, xiA3, xiA4 = self.xiA
        self.xiA = np.array([1, xiA1, 0, xiA3, 0], dtype=DOUBLE)
        self.calc_ABDE_from_lamination_parameters()


    def force_symmetric_LP(self):
        r"""Force symmetric lamination parameters

        The lamination parameters `\xi_{Bi}` are set to null
        to force a symmetric laminate.

        """
        self.xiB = np.zeros(5)
        self.calc_ABDE_from_lamination_parameters()


    def force_orthotropic(self):
        r"""Force an orthotropic laminate

        The terms
        `A_{13}`, `A_{23}`, `A_{31}`, `A_{32}`,
        `B_{13}`, `B_{23}`, `B_{31}`, `B_{32}`,
        `D_{13}`, `D_{23}`, `D_{31}`, `D_{32}` are set to zero to force an
        orthotropic laminate.

        """
        if self.offset != 0.:
            raise RuntimeError(
                    'Laminates with offset cannot be forced orthotropic!')
        self.A[0, 2] = 0.
        self.A[1, 2] = 0.
        self.A[2, 0] = 0.
        self.A[2, 1] = 0.

        self.B[0, 2] = 0.
        self.B[1, 2] = 0.
        self.B[2, 0] = 0.
        self.B[2, 1] = 0.

        self.D[0, 2] = 0.
        self.D[1, 2] = 0.
        self.D[2, 0] = 0.
        self.D[2, 1] = 0.

        self.ABD[0, 2] = 0. # A16
        self.ABD[1, 2] = 0. # A26
        self.ABD[2, 0] = 0. # A61
        self.ABD[2, 1] = 0. # A62

        self.ABD[0, 5] = 0. # B16
        self.ABD[5, 0] = 0. # B61
        self.ABD[1, 5] = 0. # B26
        self.ABD[5, 1] = 0. # B62

        self.ABD[3, 2] = 0. # B16
        self.ABD[2, 3] = 0. # B61
        self.ABD[4, 2] = 0. # B26
        self.ABD[2, 4] = 0. # B62

        self.ABD[3, 5] = 0. # D16
        self.ABD[4, 5] = 0. # D26
        self.ABD[5, 3] = 0. # D61
        self.ABD[5, 4] = 0. # D62

        self.ABDE[0, 2] = 0. # A16
        self.ABDE[1, 2] = 0. # A26
        self.ABDE[2, 0] = 0. # A61
        self.ABDE[2, 1] = 0. # A62

        self.ABDE[0, 5] = 0. # B16
        self.ABDE[5, 0] = 0. # B61
        self.ABDE[1, 5] = 0. # B26
        self.ABDE[5, 1] = 0. # B62

        self.ABDE[3, 2] = 0. # B16
        self.ABDE[2, 3] = 0. # B61
        self.ABDE[4, 2] = 0. # B26
        self.ABDE[2, 4] = 0. # B62

        self.ABDE[3, 5] = 0. # D16
        self.ABDE[4, 5] = 0. # D26
        self.ABDE[5, 3] = 0. # D61
        self.ABDE[5, 4] = 0. # D62


    def force_symmetric(self):
        """Force a symmetric laminate

        The `B` terms of the constitutive matrix are set to zero.

        """
        if self.offset != 0.:
            raise RuntimeError(
                    'Laminates with offset cannot be forced symmetric!')
        self.B = np.zeros((3,3))
        self.ABD[0:3, 3:6] = 0
        self.ABD[3:6, 0:3] = 0

        self.ABDE[0:3, 3:6] = 0
        self.ABDE[3:6, 0:3] = 0

