class Representation():
    def _repr(self, keys, skip=4*' '):
        ss="Class : {0} \n".format(self.__class__.__name__)
        nstr=np.array([key.__len__() for key in keys]).max()

        for key in keys:
            attr=getattr(self, key)
            if isinstance(attr, collections.Callable):
                ss+='{0}{1}{3} = {2}\n'.format(skip, key, str(attr()), (nstr-key.__len__())*' ')
            else:
                ss+='{0}{1}{3} = {2}\n'.format(skip, key, str(attr), (nstr-key.__len__())*' ')
        return ss class Grid():
class TrigPolBasis(Grid):
    def get_ZNl(N, fft_form=fft_form_default):
    def get_xil(N, Y, fft_form=fft_form_default):
    def get_freq(N, Y, fft_form=fft_form_default):
    def get_product(xi):
    def get_coordinates(N, Y):
    def eval_phi_k_N(self, x):
    def eval_phi_k(self, x):
    def get_nodes(self):
    def __init__(self, order, N=None, Y=None):
    def __repr__(self):
def get_inverse(A):
def enlarge(xN, M):
def decrease(xN, M):
def get_Nodd(N):
def mean_index(N, fft_form=fft_form_default):
        N = np.atleast_1d(np.array(N, dtype=np.int))
    B = np.copy(A)
    N = np.array(A.shape[2:])
    M = np.array(M, dtype=np.float)
    N = np.array(xN.shape, dtype=np.float)
    M = np.array(M, dtype=np.float)
    N = np.array(xN.shape, dtype=np.float) class ElasticTensor():
    """
    This class represents linear fourth-order tensor of elastic parameters
    for both stiffness and compliance. It also evaluates the tensor as a matrix
    using engineering Mandel's or Voigt's notation. The plane stress or strain
    is also available.
    """

    def __init__(self, bulk=None, mu=None, stiffness=True, plane=None):
        """
        Parameters
        ----------
        bulk : float
            bulk modulus
        mu : float
            shear modulus
        stiffness : boolean
            determines whether the values are evaluated for stiffness (True)
            or for its inverse compliance (False)
        plane : None, 'strain', or 'stress'
            determines the physical problem; None represents three-dimensional
            problem and 'strain' or 'stress' represent reduced problem of plane
            strain or stress respectively
        """
        if stiffness:
            self.stiffness = stiffness
            self.val_type = 'stiffness'
        else:
            self.val_type = 'compliance'

        if plane is None:
            self.dim = 3
        elif plane in ['stress', 'strain']:
            self.dim = 2
        else:
            raise ValueError("This type of plane (%s) is not supported."
                             % str(plane))

        self.plane = plane

        self.sym = self.get_sym(self.dim)
        self.bulk = float(bulk)
        self.mu = float(mu)

        # set up values to tensor and matrices
        _, volumetric, deviatoric = self.get_decomposition()
        if stiffness:
            self.val = 3.*bulk*volumetric + 2.*mu*deviatoric
        else:
            self.val = 1./(3*bulk)*volumetric + 1./(2*mu)*deviatoric

        self.mandel = self.create_mandel(self.val)
        self.voigt = self.create_voigt(self.val)

        if plane is not None:
            self.val = self.val[0:2, 0:2, 0:2, 0:2]

        if (stiffness and plane == 'strain') or \
                (not stiffness and plane == 'stress'):
            self.mandel = self.get_plane_in_engineering(self.mandel)
            self.voigt = self.get_plane_in_engineering(self.voigt)
        elif (not stiffness and plane == 'strain') or \
                (stiffness and plane == 'stress'):
            inv = np.linalg.inv
            self.mandel = inv(self.get_plane_in_engineering(inv(self.mandel)))
            self.voigt = inv(self.get_plane_in_engineering(inv(self.voigt)))
        else:
            pass

    @staticmethod
    def get_plane_in_engineering(val, ind=None):
        """
        get plane (strain or stress) of matrix in Mandel or Voigth notation

        Parameters:
        val : numpy.ndarray of shape (6,6)
            matrix to transfer to plane
        ind : list or tuple of len 2
            indicies of plane axis

        Returns:
        mat : numpy.ndarray of shape (3,3)
            plane matrix
        """
        if ind is None:
            ind = [0, 1]
        else:
            ind = list(ind)
        ind_shear = list(range(3))
        ind_shear.remove(ind[0])
        ind_shear.remove(ind[1])
        ind.append(ind_shear[0]+3)
        mat = val[ind][:, ind]
        return mat

    @staticmethod
    def get_plane_in_tensor(val, ind=None):
        assert(isinstance(val, np.ndarray))
        if ind is None:
            ind = [0, 1]

        if val.shape[:4]==4*(3,): # dimension = 4
            mat = val[ind][:,ind][:,:,ind][:,:,:,ind]
        elif val.shape[:2]==2*(3,): # dimension = 2
            mat = val[ind][:, ind]
        else:
            raise ValueError('Shape {0} of input!'.format(val.shape))
        return mat

    def __repr__(self):
        ss = "Class: %s\n" % (self.__class__.__name__)
        ss += '    stiffness = %s (%s)\n' % (self.stiffness, self.val_type)
        ss += '    dim = %d' % (self.dim)
        if self.plane is None:
            ss += '\n'
        else:
            ss += ' (plane %s)\n' % self.plane
        ss += '    bulk = %s\n' % str(self.bulk)
        ss += '    mu = %s\n' % str(self.mu)
        return ss

    @staticmethod
    def get_sym(dim):
        return dim*(dim+1)/2

    @staticmethod
    def get_decomposition():
        """
        It produces symmetrized fourth-order identity, hydrostatic, and
        deviatoric projections.

        Returns
        -------
        idsym : numpy.array of shape = (3, 3, 3, 3)
            symmetrized identity operator
        volumetric : numpy.array of shape = (3, 3, 3, 3)
            hydrostatic projection
        deviatoric : numpy.array of shape = (3, 3, 3, 3)
            deviatoric projection
        """
        ids = np.eye(3)
        volumetric = 1./3*np.einsum('ij,kl', ids, ids)
        idsym = 0.5*(np.einsum('ik,jl', ids, ids)
                     + np.einsum('il,jk', ids, ids))
        deviatoric = idsym - volumetric
        return idsym, volumetric, deviatoric

    @staticmethod
    def create_mandel(mat, ndim=None):
        """
        It transfer symmetric four-order tensor (or matrix) to
        second order tensor (or vector) using Mandel's notation.

        Parameters
        ----------
        mat : numpy.array of shape = (d, d, d, d) or (d, d) for dimension d
            fourth-order tensor of elastic parameters
        ndim : 2 or 4
            dimensionality of input

        Returns
        -------
        res : numpy.array of shape = (sym, sym) or (sym,) for sym = d*(d+1)/2
            second-order tensor of elastic parameters with Mandel's notation
        """
        dim = mat.shape[0]
        sym = int(dim*(dim+1)/2)
        if ndim is None:
            ndim = mat.ndim
        grid_shape = mat.shape[ndim:]

        if ndim == 4:
            res = np.zeros((sym, sym) + grid_shape, dtype=mat.dtype)
            if dim == 3:
                for ii, jj in np.ndindex(dim, dim):
                    kk = list(range(dim))
                    kk.remove(ii)
                    ll = list(range(dim))
                    ll.remove(jj)
                    res[ii, jj] = mat[ii, ii, jj, jj]
                    res[ii, jj+dim] = 2**.5*mat[ii, ii, ll[0], ll[1]]
                    res[jj+dim, ii] = res[ii, jj+dim]
                    res[ii+dim, jj+dim] = 2*mat[kk[0], kk[1], ll[0], ll[1]]
            if dim==2:
                res[2,2] = 2*mat[0,1,0,1]
                for ii in range(dim):
                    res[ii,2] = 2**.5*mat[ii,ii,0,1]
                    res[2,ii] = 2**.5*mat[0,1,ii,ii]
                    for jj in range(dim):
                        res[ii, jj] = mat[ii, ii, jj, jj]

        elif ndim == 2:
            res = np.zeros((sym,)  + grid_shape, dtype=mat.dtype)
            res[:dim] = np.diag(mat)
            if dim == 2:
                res[dim] = 2**.5*mat[0, 1]
            elif dim == 3:
                for ii in np.arange(sym-dim):
                    ind = list(range(sym-dim))
                    ind.remove(ii)
                    res[dim+ii] = 2**.5*mat[ind[0], ind[1]]
            else:
                raise ValueError("Incorrect dimension (%d)" % dim)
        return res

    @staticmethod
    def dispose_mandel(vec, ndim=2):
        assert(isinstance(vec, np.ndarray))
        vec = vec.squeeze()
        sym = vec.shape[0]

        def dimfun(sym):
            """
            Inverse function to dim*(dim+1)/2.
            """
            return int((-1.+(1+8*sym)**.5)/2)

        dim = dimfun(sym)
        if ndim is None:
            ndim = vec.ndim
        grid_shape = vec.shape[ndim:]

        if ndim == 2: # matrix -> tensor4
            assert(vec.shape[0] == vec.shape[1])

            mat = np.zeros(4*(dim,)+grid_shape, dtype=vec.dtype)

            if dim==3:
                for ii in np.arange(dim):
                    for jj in np.arange(dim):
                        kk = list(range(dim))
                        kk.remove(ii)
                        ll = list(range(dim))
                        ll.remove(jj)
                        mat[ii, ii, jj, jj] = vec[ii, jj]
                        mat[ii, ii, ll[0], ll[1]] = vec[ii, jj+dim] / 2**.5
                        mat[ii, ii, ll[1], ll[0]] = vec[ii, jj+dim] / 2**.5
                        mat[ll[0], ll[1], ii, ii] = mat[ii, ii, ll[0], ll[1]]
                        mat[ll[1], ll[0], ii, ii] = mat[ii, ii, ll[0], ll[1]]
                        mat[kk[0], kk[1], ll[0], ll[1]] = vec[ii+dim, jj+dim] / 2.
                        mat[kk[1], kk[0], ll[0], ll[1]] = vec[ii+dim, jj+dim] / 2.
                        mat[kk[0], kk[1], ll[1], ll[0]] = vec[ii+dim, jj+dim] / 2.
                        mat[kk[1], kk[0], ll[1], ll[0]] = vec[ii+dim, jj+dim] / 2.
            elif dim==2:
                mat[0,1,0,1] = vec[2,2]/2.
                mat[0,1,1,0] = vec[2,2]/2.
                mat[1,0,0,1] = vec[2,2]/2.
                mat[1,0,1,0] = vec[2,2]/2.
                for ii in range(dim):
                    mat[ii,ii,0,1] = vec[ii,2]/2**.5
                    mat[ii,ii,1,0] = vec[ii,2]/2**.5
                    mat[0,1,ii,ii] = vec[2,ii]/2**.5
                    mat[1,0,ii,ii] = vec[2,ii]/2**.5
                    for jj in range(dim):
                        mat[ii,ii,jj,jj] = vec[ii,jj]
                        mat[ii,ii,jj,jj] = vec[ii,jj]

        elif ndim == 1: # vector -> matrix
            mat = np.diag(vec[:dim])

            if dim == 2:
                mat[0, 1] = vec[-1]/2**0.5
                mat[1, 0] = vec[-1]/2**0.5
            elif dim == 3:
                for ii in np.arange(sym-dim):
                    ind = list(range(sym-dim))
                    ind.remove(ii)
                    mat[ind[0], ind[1]] = vec[dim+ii]/2.**.5
                    mat[ind[1], ind[0]] = vec[dim+ii]/2.**.5
            else:
                raise ValueError("Incorrect dimension (%d)" % dim)

        return mat

    @staticmethod
    def create_voigt(mat, valtype='strain'):
        """
        It transfer symmetric four-order tensor to second order tensor
        using Voigt's notation.

        Parameters
        ----------
        mat : numpy.array of shape = (3, 3, 3, 3)
            fourth-order tensor of elastic parameters
        valtype : one of 'strain' or 'stress'
            this distinguish a engineering notation for strain and stress

        Returns
        -------
        vec : numpy.array of shape = (6, 6)
            second-order tensor of elastic parameters with Voigt's notation
        """
        dim = mat.shape[0]
        sym = int(dim*(dim+1)/2)
        if mat.ndim == 4:
            vec = np.zeros([sym, sym], dtype=mat.dtype)
            for ii in np.arange(dim):
                for jj in np.arange(dim):
                    kk = list(range(dim))
                    kk.remove(ii)
                    ll = list(range(dim))
                    ll.remove(jj)
                    vec[ii, jj] = mat[ii, ii, jj, jj]
                    vec[ii, jj+dim] = mat[ii, ii, ll[0], ll[1]]
                    vec[jj+dim, ii] = vec[ii, jj+dim]
                    vec[ii+dim, jj+dim] = mat[kk[0], kk[1], ll[0], ll[1]]

        elif mat.ndim == 2:
            vec = np.zeros(sym, dtype=mat.dtype)
            vec[:dim] = np.diag(mat)
            if valtype == 'strain':
                coef = 2.
            elif valtype == 'stress':
                coef = 1.
            else:
                msg = "Parameter valtype (%s) should be one of 'strain' or\
                     'stress'." % (str(valtype),)
                raise ValueError(msg)

            if dim == 2:
                vec[dim] = coef*mat[0, 1]
            elif dim == 3:
                for ii in np.arange(sym-dim):
                    ind = list(range(sym-dim))
                    ind.remove(ii)
                    vec[dim+ii] = coef*mat[ind[0], ind[1]]
            else:
                raise ValueError("Incorrect dimension (%d)" % dim)
        return vec def enlarge(xN, M):
    """
    Enlarge an array of Fourier coefficients by zeros.

    Parameters
    ----------
    xN : numpy.ndarray of shape = N
        input array that is to be enlarged

    Returns
    -------
    xM : numpy.ndarray of shape = M
        output array that is enlarged
    M : array like
        number of grid points
    """
    M = np.array(M, dtype=np.float)
    N = np.array(xN.shape, dtype=np.float)
    if np.allclose(M, N):
        return xN

    ibeg = np.ceil((M-N)/2).astype(dtype=np.int)
    iend = np.ceil((M+N)/2).astype(dtype=np.int)

    slc=[slice(ibeg[i],iend[i],1) for i in range(N.size)]
    xM = np.zeros(M.astype(dtype=np.int), dtype=xN.dtype)
    xM[tuple(slc)]=xN
    return xM def decrease(xN, M):
    """
    Decreases an array of Fourier coefficients by omitting the highest
    frequencies.

    Parameters
    ----------
    xN : numpy.ndarray of shape = N
        input array that is to be enlarged

    Returns
    -------
    xM : numpy.ndarray of shape = M
        output array that is enlarged
    M : array like
        number of grid points
    """
    M = np.array(M, dtype=np.float)
    N = np.array(xN.shape, dtype=np.float)
    ibeg = np.fix((N-M+(M % 2))/2).astype(dtype=np.int)
    iend = np.fix((N+M+(M % 2))/2).astype(dtype=np.int)

    slc=[slice(ibeg[i],iend[i],1) for i in range(N.size)]
    return xN[tuple(slc)] def get_inverse(A):
    """
    It calculates the inverse of conductivity coefficients at grid points,
    i.e. of matrix A_GaNi

    Parameters
    ----------
    A : numpy.ndarray

    Returns
    -------
    invA : numpy.ndarray
    """
    B = np.copy(A)
    N = np.array(A.shape[2:])
    d = A.shape[0]
    if A.shape[0] != A.shape[1]:
        raise NotImplementedError("Non-square matrix!")

    invA = np.eye(d).tolist()
    for m in np.arange(d):
        Bdiag = np.copy(B[m][m])
        B[m][m] = np.ones(N)
        for n in np.arange(m+1, d):
            B[m][n] = B[m][n]/Bdiag
        for n in np.arange(d):
            invA[m][n] = invA[m][n]/Bdiag
        for k in np.arange(m+1, d):
            Bnull = np.copy(B[k][m])
            for l in np.arange(d):
                B[k][l] = B[k][l] - B[m][l]*Bnull
                invA[k][l] = invA[k][l] - invA[m][l]*Bnull
    for m in np.arange(d-1, -1, -1):
        for k in np.arange(m-1, -1, -1):
            Bnull = np.copy(B[k][m])
            for l in np.arange(d):
                B[k][l] = B[k][l] - B[m][l]*Bnull
                invA[k][l] = invA[k][l] - invA[m][l]*Bnull
    invA = np.array(invA)
    return invA class Grid():
    @staticmethod
    def get_ZNl(N, fft_form=fft_form_default):
        r"""
        it produces index set ZNl=\underline{\set{Z}}^d_N :
        ZNl[i][j]\in\set{Z} : -N[i]/2 <= ZNl[i] < N[i]/2
        """
        ZNl = []
        N = np.atleast_1d(np.array(N, dtype=np.int))
        for m in range(N.size):
            ZNl.append(np.arange(np.fix(-N[m]/2.), np.fix(N[m]/2.+0.5),
                                 dtype=np.int))
        if fft_form in ['r',0]:
            return [np.fft.ifftshift(val) for val in ZNl]
        else:
            return ZNl

    @staticmethod
    def get_xil(N, Y, fft_form=fft_form_default):
        """
        it produces discrete frequencies of Fourier series
        xil[i] = ZNl[i]/Y[i]
        """
        xil = []
        for m in np.arange(np.size(N)):
            xil.append(np.arange(np.fix(-N[m]/2.), np.fix(N[m]/2.+0.5))/Y[m])
        if fft_form in ['r',]:
            xil=[np.fft.ifftshift(xi) for xi in xil]
            xil[-1] = xil[-1][:int(np.fix(N[-1]/2)+1)]
        elif fft_form in [0]:
            xil = [np.fft.ifftshift(xi) for xi in xil]
        return xil

    @staticmethod
    def get_freq(N, Y, fft_form=fft_form_default):
        return Grid.get_xil(N, Y, fft_form=fft_form)

    @staticmethod
    def get_product(xi):
        xis = np.atleast_2d(xi[0])
        for ii in range(1, len(xi)):
            xis_new = np.tile(xi[ii], xis.shape[1])
            xis_old = np.repeat(xis, xi[ii].size, axis=1)
            xis = np.vstack([xis_old, xis_new])
        return xis

    @staticmethod
    def get_coordinates(N, Y):
        """
        It produces coordinates of the set of nodal points
        Coord[i][j] = x_N^{(i,j)}
        """
        d = np.size(N)
        ZNl = Grid.get_ZNl(N, fft_form='c')
        coord = np.zeros(np.hstack([d, N]))
        for ii in np.arange(d):
            x = Y[ii]*ZNl[ii]/N[ii]
            Nshape = np.ones(d, dtype=np.int)
            Nshape[ii] = N[ii]
            Nrep = np.copy(N)
            Nrep[ii] = 1
            coord[ii] = np.tile(np.reshape(x, Nshape), Nrep)
        return coord def fftn(x, N): # normalised FFT
    return 1./np.prod(N)*fft.fftn(x, N) def ifftn(x, N): # normalised FFT
    return fft.ifftn(x, N).real*np.prod(N) def fftnc(x, N):
    """
    Fourier centered n-dimensional FFT algorithm
    """
    ax=tuple(np.setdiff1d(list(range(x.ndim)), list(range(x.ndim-N.__len__())), assume_unique=True))
    return 1./np.prod(N)*fft.fftshift(fft.fftn(x, N), ax) def icfftn(Fx, N):
    """
    Fourier centered n-dimensional inverse FFT algorithm
    """
    ax=tuple(np.setdiff1d(list(range(Fx.ndim)), list(range(Fx.ndim-N.__len__())), assume_unique=True))
    return fft.ifftn(fft.ifftshift(Fx, ax), N).real*np.prod(N) def rfftn(x, N): # real-valued FFT
    return fft.rfftn(x, N) def irfftn(x, N): # real-valued FFT
    return fft.irfftn(x, N)
import numpy as np
import itertools
        import matplotlib.pyplot as plt
from ffthompy.general.base import Representation
from ffthompy.trigpol import mean_index, fft_form_default, get_Nodd
from ffthompy.mechanics.matcoef import ElasticTensor
from ffthompy.trigpol import enlarge, decrease, get_inverse, Grid
from ffthompy.tensors.fft import fftn, ifftn, fftnc, icfftn, rfftn, irfftn
from copy import copy
        from mpl_toolkits.mplot3d import axes3d
        from ffthompy.trigpol import Grid
            from matplotlib import cm
                self.val[di]+=mean[di]
        return self

    def __eq__(self, Y, full=True, tol=1e-13):
        """
        Check the equality with other objects comparable to trig. polynomials.
        """
        X=self
        _bool=False
        res=np.inf
        if (isinstance(X, Tensor) and X.fft_form==Y.fft_form and
                X.val.squeeze().shape==Y.val.squeeze().shape and X.Fourier==Y.Fourier):
            res=np.linalg.norm(X.val.squeeze()-Y.val.squeeze())
            if res<tol:
                _bool=True
        if full:
            return _bool, res
        else:
            return _bool

    def set_shape(self):
        shape_size=self.val.ndim-self.N.size
        self.shape=np.array(self.val.shape[:shape_size])
        return self.shape

    def transpose(self):
        if self.order==2:
            val=np.einsum('ij...->ji...', self.val)
        elif self.order==4:
            val=np.einsum('ijkl...->klij...', self.val)
        else:
            raise NotImplementedError()
        return self.copy(name=self.name[:10]+'.T', val=val)

    def transpose_left(self):
        res=self.empty_like(name=self.name[:10]+'.T')
        assert(self.order==4)
        res.val=np.einsum('ijkl...->jikl...', self.val)
        return res

    def transpose_right(self):
        res=self.empty_like(name=self.name[:10]+'.T')
        assert(self.order==4)
        res.val=np.einsum('ijkl...->ijlk...', self.val)
        return res

    def identity(self):
        self.val[:]=0.
        assert(self.order % 2 == 0)
        for ii in itertools.product(*tuple([list(range(n)) for n in self.shape[:int(self.order/2)]])):
            self.val[ii+ii]=1.

    def vec(self):
        """
        Returns one-dimensional vector (column) version of trigonometric
        polynomial.
        """
        return np.matrix(self.val.ravel()).transpose()

    def zeros_like(self, name=None):
        if name is None:
            name='zeros({})'.format(self.name[:10])
        return self.copy(name=name, val=np.zeros_like(self.val))

    def empty_like(self, name=None):
        if name is None:
            name='empty({})'.format(self.name[:10])
        return self.copy(name=name, val=np.empty_like(self.val))

    def calc_eigs(self, sort=True, symmetric=False, mandel=False):
        if symmetric:
            eigfun=np.linalg.eigvalsh
        else:
            eigfun=np.linalg.eigvals

        if self.order==2:
            eigs=np.zeros(self.N+(self.shape[0],))
            for ind in np.ndindex(self.N):
                mat=self.val[(slice(None), slice(None))+ind]
                eigs[ind]=eigfun(mat)

        elif self.order==4:
            if mandel:
                matrixfun=lambda x: ElasticTensor.create_mandel(x)
                d=self.shape[2]
                eigdim=d*(d+1)/2
            else:
                eigdim=self.shape[2]*self.shape[3]
                matrixfun=lambda x: np.reshape(x, 2*(eigdim,))

            eigs=np.zeros(self.N+(eigdim,))
            val=np.copy(self.val)
            for ii in range(self.dim):
                val=np.rollaxis(val, self.val.ndim-self.dim+ii, ii)

            for ind in np.ndindex(*self.N):
                eigs[ind]=eigfun(matrixfun(val[ind]))

        eigs=np.rollaxis(np.array(eigs), -1, 0)
        if sort:
            eigs=np.sort(eigs, axis=0)
        return eigs

    @property
    def axes(self): # axes for Fourier transform
        return tuple(range(self.order, self.order+self.dim))

    def fourier(self, Fourier=None, copy=False):
        assert(self.origin==0)

        if self.Fourier==Fourier:
            if copy:
                return self.copy()
            else:
                return self

        if copy:
            if self.Fourier:
                return self.copy(val=self.ifftn(self.val, self.N), Fourier=not self.Fourier)
            else:

                return self.copy(val=self.fftn(self.val, self.N), Fourier=not self.Fourier)