from __future__ import division, absolute_import

import platform
import gc
import pickle
from multiprocessing import cpu_count

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs, eigsh
from scipy.linalg import eig
from numpy import linspace, deg2rad

import compmech.composite.laminate as laminate
from compmech.analysis import Analysis
from compmech.logger import msg, warn
from compmech.constants import DOUBLE
from compmech.sparse import (make_symmetric, remove_null_cols,
        make_skew_symmetric)

from . import modelDB


def load(name):
    if '.Panel' in name:
        return pickle.load(open(name, 'rb'))
    else:
        return pickle.load(open(name + '.Panel', 'rb'))


def check_c(c, size):
    if not isinstance(c, np.ndarray):
        raise TypeError('"c" must be a NumPy ndarray object')
    if c.ndim != 1:
        raise ValueError('"c" must be a 1-D ndarray object')
    if c.shape[0] != size:
        raise ValueError('"c" must have the same size as the global stiffness matrix')


class Panel(object):
    r"""General Panel class

    It works for both flat plates, cylindrical and conical panels. The right
    model is selected according to parameters ``r`` (radius) and ``alphadeg``
    (semi-vertex angle).

    The approximation functions for the displacement fields are built using
    :ref:`Bardell's functions <theory_func_bardell>`.

    Parameters
    ----------
    a : float, optional
        Length (along the `x` coordinate).
    b : float, optional
        Width (along the `y` coordinate).
    r : float, optional
        Radius for cylindrical panels.
    alphadeg : float, optional
        Semi-vertex angle for conical panels.
    stack : list or tuple, optional
        A sequence representing the angles for each ply.
    plyt : float, optional
        Ply thickness.
    laminaprop : list or tuple, optional
        Orthotropic lamina properties: `E_1, E_2, \nu_{12}, G_{12}, G_{13}, G_{23}`.
    mu : float, optional
        Material density.
    m, n : int, optional
        Number of terms for the approximation functions along `x` and `y`,
        respectively.
    offset : float, optional
        Laminate offset about panel mid-surface. The offset is measured along
        the normal (`z`) axis.
    y1, y2 : float, optional
        Define the lower and upper limit along `y` for panel with incomplete
        domains.

    """
    def __init__(self, a=None, b=None, r=None, alphadeg=None,
            stack=None, plyt=None, laminaprop=None, mu=None, m=11, n=11,
            offset=0., y1=None, y2=None, **kwargs):
        self.a = a
        self.b = b
        self.y1 = y1
        self.y2 = y2
        self.r = r
        self.alphadeg = alphadeg
        self.stack = stack
        self.plyt = plyt
        self.laminaprop = laminaprop
        self.offset = offset
        # assembly
        self.group = None
        self.x0 = None
        self.y0 = None

        self.name = 'panel'
        self.bay = None

        # model
        self.model = None
        self.fsdt_shear_correction = 5/6. # in case of First-order Shear Deformation Theory

        # approximation series
        self.m = m
        self.n = n

        # numerical integration
        self.nx = m
        self.ny = n
        self.ni_num_cores = cpu_count()//2
        self.ni_method = 'trapz2d'
        self.c0 = None

        # loads
        self.Nxx = None
        self.Nyy = None
        self.Nxy = None
        self.Nxx_cte = None
        self.Nyy_cte = None
        self.Nxy_cte = None
        self.forces = []
        self.forces_inc = []

        # boundary conditions

        # displacement at 4 edges is zero
        # free to rotate at 4 edges (simply supported by default)
        self.u1tx = 0.
        self.u1rx = 0.
        self.u2tx = 0.
        self.u2rx = 0.
        self.v1tx = 0.
        self.v1rx = 0.
        self.v2tx = 0.
        self.v2rx = 0.
        self.w1tx = 0.
        self.w1rx = 1.
        self.w2tx = 0.
        self.w2rx = 1.

        self.u1ty = 0.
        self.u1ry = 0.
        self.u2ty = 0.
        self.u2ry = 0.
        self.v1ty = 0.
        self.v1ry = 0.
        self.v2ty = 0.
        self.v2ry = 0.
        self.w1ty = 0.
        self.w1ry = 1.
        self.w2ty = 0.
        self.w2ry = 1.

        # material
        self.mu = mu
        self.plyts = None
        self.laminaprops = None

        # aeroelastic parameters
        self.flow = 'x'
        self.beta = None
        self.gamma = None
        self.aeromu = None
        self.rho_air = None
        self.speed_sound = None
        self.Mach = None
        self.V = None

        # constitutive law
        self.F = None
        self.force_orthotropic_laminate = False

        # eigenvalue analysis
        self.num_eigvalues = 5
        self.num_eigvalues_print = 5

        # output queries
        self.out_num_cores = cpu_count()

        # analysis
        self.analysis = Analysis(self.calc_fext, self.calc_k0, self.calc_fint, self.calc_kT)

        # outputs
        self.increments = None
        self.eigvecs = None
        self.eigvals = None

        for k, v in kwargs.items():
            setattr(self, k, v)

        self._clear_matrices()


    def _clear_matrices(self):
        self.k0 = None
        self.kG0 = None
        self.kM = None
        self.kA = None
        self.cA = None
        self.lam = None
        self.u = None
        self.v = None
        self.w = None
        self.phix = None
        self.phiy = None
        self.Xs = None
        self.Ys = None

        #NOTE memory cleanup
        gc.collect()


    def _rebuild(self):
        if self.model is None:
            if self.r is None and self.alphadeg is None:
                self.model = 'plate_clt_donnell_bardell'
            elif self.r is not None and self.alphadeg is None:
                self.model = 'cpanel_clt_donnell_bardell'
            elif self.r is not None and self.alphadeg is not None:
                self.model = 'kpanel_clt_donnell_bardell'

        valid_models = sorted(modelDB.db.keys())

        if not self.model in valid_models:
            raise ValueError('ERROR - valid models are:\n    ' +
                     '\n    '.join(valid_models))

        if not self.stack:
            raise ValueError('stack must be defined')

        if not self.laminaprops:
            if not self.laminaprop:
                raise ValueError('laminaprop must be defined')
            self.laminaprops = [self.laminaprop for i in self.stack]

        if not self.plyts:
            if self.plyt is None:
                raise ValueError('plyt must be defined')
            self.plyts = [self.plyt for i in self.stack]


    def get_size(self):
        r"""Calculates the size of the stiffness matrices

        The size of the stiffness matrices can be interpreted as the number of
        rows or columns, recalling that this will be the size of the Ritz
        constants' vector `\{c\}`, the internal force vector `\{F_{int}\}` and
        the external force vector `\{F_{ext}\}`.

        Returns
        -------
        size : int
            The size of the stiffness matrices.

        """
        num = modelDB.db[self.model]['num']
        self.size = num*self.m*self.n
        return self.size


    def _default_field(self, xs, ys, gridx, gridy):
        if xs is None or ys is None:
            xs = linspace(0, self.a, gridx)
            ys = linspace(0, self.b, gridy)
            xs, ys = np.meshgrid(xs, ys, copy=True)
        xs = np.atleast_1d(np.array(xs, dtype=DOUBLE))
        ys = np.atleast_1d(np.array(ys, dtype=DOUBLE))
        xshape = xs.shape
        tshape = ys.shape
        if xshape != tshape:
            raise ValueError('Arrays xs and ys must have the same shape')
        self.Xs = xs
        self.Ys = ys
        xs = np.ascontiguousarray(xs.ravel(), dtype=DOUBLE)
        ys = np.ascontiguousarray(ys.ravel(), dtype=DOUBLE)

        return xs, ys, xshape, tshape


    def _get_lam_F(self):
        if self.lam is None:
            raise RuntimeError('lam object is None!')
        if 'clt' in self.model:
            F = self.lam.ABD

        elif 'fsdt' in self.model:
            F = self.lam.ABDE
            F[6:, 6:] *= self.fsdt_shear_correction

        if self.force_orthotropic_laminate:
            msg('', silent=silent)
            msg('Forcing orthotropic laminate...', level=2, silent=silent)
            F[0, 2] = 0. # A16
            F[1, 2] = 0. # A26
            F[2, 0] = 0. # A61
            F[2, 1] = 0. # A62

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

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

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

            if F.shape[0] == 8:
                F[6, 7] = 0. # A45
                F[7, 6] = 0. # A54

        return F


    def calc_k0(self, size=None, row0=0, col0=0, silent=False, finalize=True,
            c=None, nx=None, ny=None, Fnxny=None, inc=None, NLgeom=False):
        """Calculate the constitutive stiffness matrix

        If ``c`` is not given it calculates the linear constitutive stiffness
        matrix, otherwise the large displacement linear constitutive stiffness
        matrix is calculated. When using ``c``the size of ``c`` must be the
        same as ``size``.

        In assemblies of semi-analytical models the sparse matrices that are
        calculated may have the ``size`` of the assembled global model, and the
        current constitutive matrix being calculated starts at position
        ``row0`` and ``col0``.

        Parameters
        ----------
        size : int
            The size of the calculated sparse matrices.
        row0, col0: int or None, optional
            Offset to populate the output sparse matrix (useful when
            assemblying panels).
        silent : bool, optional
            A boolean to tell whether the log messages should be printed.
        finalize : bool, optional
            Asserts validity of output data and makes the output matrix
            symmetric, should be ``False`` when assemblying.
        c : array-like or None, optional
            This must be the result of a static analysis, used to compute the
            non-linear term based on the actual displacement field.
        nx, ny : int or None, optional
            Number of integration points along `x` and `y`, respectively, for
            the Legendre-Gauss quadrature rule applied in the numerical
            integration. Only used when ``c`` is given.
        Fnxny : 4-D array-like or None, optional
            The constitutive relations for the laminate at each integration
            point. Must be a 4-D array of shape ``(nx, ny, 6, 6)`` when using
            classical laminated plate theory models.
        NLgeom : bool, optional
            Flag to indicate if geometrically non-linearities should be
            considered.

        """
        self._rebuild()
        if size is None:
            size = self.get_size()
        if c is None:
            msg('Calculating k0... ', level=2, silent=silent)
        else:
            check_c(c, size)
            msg('Calculating kL... ', level=2, silent=silent)

        matrices = modelDB.db[self.model]['matrices']

        alphadeg = self.alphadeg if self.alphadeg is not None else 0.
        self.alpharad = deg2rad(alphadeg)
        self.r = self.r if self.r is not None else 0.

        if self.stack is not None:
            lam = laminate.read_stack(self.stack, plyts=self.plyts,
                                      laminaprops=self.laminaprops,
                                      offset=self.offset)
            self.lam = lam
            self.F = self._get_lam_F()

        if self.y1 is not None and self.y2 is not None:
            if c is not None or Fnxny is not None:
                raise NotImplementedError(
                'Partial domain from y1 to y2 not implemented for kL')
            k0 = matrices.fk0y1y2(self.y1, self.y2, self, size, row0, col0)
        else:
            if c is None and Fnxny is None:
                k0 = matrices.fk0(self, size, row0, col0)
            else:
                matrices_num = modelDB.db[self.model]['matrices_num']
                nx = self.nx if nx is None else nx
                ny = self.ny if ny is None else ny
                #NOTE the consistence checks for Fnxny are done within the .pyx
                #     files
                Fnxny = self.F if Fnxny is None else Fnxny
                if c is None:
                    # Empty c if the interest is only on the heterogeneous
                    # laminate properties
                    c = np.zeros(self.size, dtype=DOUBLE)
                c = np.ascontiguousarray(c, dtype=DOUBLE)
                k0 = matrices_num.fkL_num(c, Fnxny, self,
                         size, row0, col0, nx, ny, NLgeom=int(NLgeom))

        #TODO allow constant stress state to be obtained using static results,
        #     which would require just passing a 'c_cte' used to calculate kG0
        Nxx_cte = self.Nxx_cte if self.Nxx_cte is not None else 0.
        Nyy_cte = self.Nyy_cte if self.Nyy_cte is not None else 0.
        Nxy_cte = self.Nxy_cte if self.Nxy_cte is not None else 0.

        if Nxx_cte != 0. or Nyy_cte != 0. or Nxy_cte != 0.:
            if self.y1 is not None and self.y2 is not None:
                k0 += matrices.fkG0y1y2(self.y1, self.y2, Nxx_cte, Nyy_cte, Nxy_cte,
                           self, size, row0, col0)
            else:
                k0 += matrices.fkG0(Nxx_cte, Nyy_cte, Nxy_cte,
                           self, size, row0, col0)

        if finalize:
            # performing final checks for k0 and making symmetry
            assert np.any(np.isnan(k0.data)) == False
            assert np.any(np.isinf(k0.data)) == False
            k0 = csr_matrix(make_symmetric(k0))

        self.k0 = k0

        #NOTE forcing Python garbage collector to clean the memory
        #     it DOES make a difference! There is a memory leak not
        #     identified, probably in the csr_matrix process
        gc.collect()

        msg('finished!', level=2, silent=silent)

        return k0


    def calc_kG0(self, size=None, row0=0, col0=0, silent=False, finalize=True,
            c=None, nx=None, ny=None, Fnxny=None, NLgeom=False):
        """Calculate the linear geometric stiffness matrix

        See :meth:`.Panel.calc_k0` for details on each parameter.

        """
        self._rebuild()
        if size is None:
            size = self.get_size()
        if c is None:
            msg('Calculating kG0... ', level=2, silent=silent)
            matrices = modelDB.db[self.model]['matrices']
        else:
            check_c(c, size)
            msg('Calculating kG... ', level=2, silent=silent)
            matrices = modelDB.db[self.model]['matrices_num']

        y1 = self.y1
        y2 = self.y2
        alphadeg = self.alphadeg if self.alphadeg is not None else 0.
        self.alpharad = deg2rad(alphadeg)
        self.r = self.r if self.r is not None else 0.

        Nxx = self.Nxx if self.Nxx is not None else 0.
        Nyy = self.Nyy if self.Nyy is not None else 0.
        Nxy = self.Nxy if self.Nxy is not None else 0.

        if c is None:
            if y1 is not None and y2 is not None:
                kG0 = matrices.fkG0y1y2(y1, y2, Nxx, Nyy, Nxy, self,
                           size, row0, col0)
            else:
                kG0 = matrices.fkG0(Nxx, Nyy, Nxy, self, size, row0, col0)
        else:
            if y1 is not None or y2 is not None:
                raise NotImplementedError('Only y1=0, y2=b is implemented!')
            c = np.ascontiguousarray(c, dtype=DOUBLE)
            nx = self.nx if nx is None else nx
            ny = self.ny if ny is None else ny
            if Fnxny is None:
                Fnxny = self._get_lam_F()
            kG0 = matrices.fkG_num(c, Fnxny, self,
                       size, row0, col0, nx, ny, NLgeom=int(NLgeom))

        if finalize:
            assert np.any(np.isnan(kG0.data)) == False
            assert np.any(np.isinf(kG0.data)) == False
            kG0 = csr_matrix(make_symmetric(kG0))

        self.kG0 = kG0

        #NOTE memory cleanup
        gc.collect()

        msg('finished!', level=2, silent=silent)

        return kG0


    def calc_kT(self, size=None, row0=0, col0=0, silent=False, finalize=True,
            c=None, nx=None, ny=None, Fnxny=None, inc=None):
        kL = self.calc_k0(size=size, row0=row0, col0=col0, silent=silent, finalize=finalize,
            c=c, nx=nx, ny=ny, Fnxny=Fnxny, inc=inc, NLgeom=True)
        kG = self.calc_kG0(size=size, row0=row0, col0=col0, silent=silent, finalize=finalize,
            c=c, nx=nx, ny=ny, Fnxny=Fnxny, NLgeom=True)
        kT = kL + kG
        self.kT = kT

        return kT


    def calc_kM(self, size=None, row0=0, col0=0, silent=False, finalize=True):
        """Calculate the mass matrix
        """
        msg('Calculating kM... ', level=2, silent=silent)

        matrices = modelDB.db[self.model]['matrices']

        y1 = self.y1
        y2 = self.y2
        alphadeg = self.alphadeg if self.alphadeg is not None else 0.
        self.alpharad = deg2rad(alphadeg)
        self.r = self.r if self.r is not None else 0.

        if size is None:
            size = self.get_size()

        #TODO allow a distribution of mu instead of constant value, at least allow a mu for each ply
        if self.mu is None:
            raise ValueError('Attribute "mu" (density) must be defined')


        h = sum(self.plyts)
        if y1 is not None and y2 is not None:
            kM = matrices.fkMy1y2(y1, y2, self.offset, self, size, row0, col0)
        else:
            kM = matrices.fkM(self.offset, self, size, row0, col0)

        if finalize:
            assert np.any(np.isnan(kM.data)) == False
            assert np.any(np.isinf(kM.data)) == False
            kM = csr_matrix(make_symmetric(kM))

        self.kM = kM

        #NOTE memory cleanup
        gc.collect()

        msg('finished!', level=2, silent=silent)

        return kM


    def calc_kA(self, size=None, row0=0, col0=0, silent=False, finalize=True):
        """Calculate the aerodynamic matrix using the linear piston theory
        """
        msg('Calculating kA... ', level=2, silent=silent)

        if 'kpanel' in self.model:
            raise NotImplementedError('Conical panels not supported')

        matrices = modelDB.db[self.model]['matrices']

        if size is None:
            size = self.get_size()

        self.r = self.r if self.r is not None else 0.

        if self.beta is None:
            if self.Mach is None:
                raise ValueError('Mach number cannot be a NoneValue')
            elif self.Mach < 1:
                raise ValueError('Mach number must be >= 1')
            elif self.Mach == 1:
                self.Mach = 1.0001
            Mach = self.Mach
            beta = self.rho_air * self.V**2 / (Mach**2 - 1)**0.5
            if self.r != 0.:
                gamma = beta*1./(2.*self.r*(Mach**2 - 1)**0.5)
            else:
                gamma = 0.
            ainf = self.speed_sound
            aeromu = beta/(Mach*ainf)*(Mach**2 - 2)/(Mach**2 - 1)
        else:
            beta = self.beta
            gamma = self.gamma if self.gamma is not None else 0.
            aeromu = self.aeromu if self.aeromu is not None else 0.

        if self.flow.lower() == 'x':
            kA = matrices.fkAx(beta, gamma, self, size, row0, col0)
        elif self.flow.lower() == 'y':
            kA = mod.fkAy(beta, self, size, row0, col0)
        else:
            raise ValueError('Invalid flow value, must be x or y')

        if finalize:
            assert np.any(np.isnan(kA.data)) == False
            assert np.any(np.isinf(kA.data)) == False
            kA = csr_matrix(make_skew_symmetric(kA))

        self.kA = kA

        #NOTE memory cleanup
        gc.collect()

        msg('finished!', level=2, silent=silent)

        return kA


    def calc_cA(self, aeromu, silent=False, finalize=True):
        """Calculate the aerodynamic damping matrix using the piston theory
        """
        msg('Calculating cA... ', level=2, silent=silent)

        cA = matrices.fcA(aeromu, self, self.size, 0, 0)
        cA = cA*(0+1j)

        if finalize:
            assert np.any(np.isnan(cA.data)) == False
            assert np.any(np.isinf(cA.data)) == False
            cA = csr_matrix(make_symmetric(cA))

        self.cA = cA

        #NOTE memory cleanup
        gc.collect()

        msg('finished!', level=2, silent=silent)


    def lb(self, tol=0, sparse_solver=True, calc_kA=False, silent=False,
           nx=10, ny=10, c=None, ckL=None, Fnxny=None):
        """Performs a linear buckling analysis

        .. note:: This will be deprecated soon, use
                  :func:`.compmech.analysis.lb`.

        The following parameters will affect the linear buckling analysis:

        =======================    =====================================
        parameter                  description
        =======================    =====================================
        ``num_eigvalues``        Number of eigenvalues to be extracted
        ``num_eigvalues_print``    Number of eigenvalues to print after
                                   the analysis is completed
        =======================    =====================================

        Parameters
        ----------
        tol : float, optional
            A float tolerance passsed to the eigenvalue solver.
        sparse_solver : bool, optional
            Tells if solver :func:`scipy.linalg.eigh` or
            :func:`scipy.sparse.linalg.eigsh` should be used.
        calc_kA : bool, optional
            If the Aerodynamic matrix should be considered.
        silent : bool, optional
            A boolean to tell whether the log messages should be printed.
        c : array-like, optional
            A set of Ritz constants that will be use to compute KG.
        ckL : array-like, optional
            A set of Ritz constants that will be use to compute KL.
        nx and ny : int or None, optional
            Number of integration points along `x` and `y`, respectively, for
            the Legendre-Gauss quadrature rule applied in the numerical
            integration.
        Fnxny : 4-D array-like or None, optional
            The constitutive relations for the laminate at each integration
            point. Must be a 4-D array of shape ``(nx, ny, 6, 6)`` when using
            classical laminated plate theory models.

        Notes
        -----
        The extracted eigenvalues are stored in the ``eigvals`` parameter
        of the ``Panel`` object and the `i^{th}` eigenvector in the
        ``eigvecs[:, i-1]`` parameter.

        """
        msg('Running linear buckling analysis...', silent=silent)

        msg('Eigenvalue solver... ', level=2, silent=silent)

        nx = self.nx if nx is None else nx
        ny = self.ny if ny is None else ny
        if ckL is None:
            self.calc_k0(silent=silent)
        else:
            self.calc_k0(silent=silent, c=ckL, nx=nx, ny=ny, Fnxny=Fnxny)
        self.calc_kG0(silent=silent, c=c, nx=nx, ny=ny, Fnxny=Fnxny)

        if calc_kA:
            raise NotImplementedError('kA requires non-Hermitian eigen solver')
            self.calc_kA(silent=silent)
            kA = self.kA
        else:
            kA = self.k0*0

        M = self.k0 + kA
        A = self.kG0

        if sparse_solver:
            mode = 'cayley'
            try:
                msg('eigsh() solver...', level=3, silent=silent)
                eigvals, eigvecs = eigsh(A=A, k=self.num_eigvalues,
                        which='SM', M=M, tol=tol, sigma=1., mode=mode)
                msg('finished!', level=3, silent=silent)
            except Exception as e:
                warn(str(e), level=4, silent=silent)
                msg('aborted!', level=3, silent=silent)
                sizebkp = A.shape[0]
                M, A, used_cols = remove_null_cols(M, A, silent=silent)
                msg('eigsh() solver...', level=3, silent=silent)
                eigvals, peigvecs = eigsh(A=A, k=self.num_eigvalues,
                        which='SM', M=M, tol=tol, sigma=1., mode=mode)
                msg('finished!', level=3, silent=silent)
                eigvecs = np.zeros((sizebkp, self.num_eigvalues),
                                   dtype=DOUBLE)
                eigvecs[used_cols, :] = peigvecs

        else:
            from scipy.linalg import eigh

            size = A.shape[0]
            M, A, used_cols = remove_null_cols(M, A, silent=silent)
            M = M.toarray()
            A = A.toarray()
            msg('eigh() solver...', level=3, silent=silent)
            eigvals, peigvecs = eigh(a=A, b=M)
            msg('finished!', level=3, silent=silent)
            eigvecs = np.zeros((size, self.num_eigvalues), dtype=DOUBLE)
            eigvecs[used_cols, :] = peigvecs[:, :self.num_eigvalues]

        eigvals = -1./eigvals

        self.eigvals = eigvals
        self.eigvecs = eigvecs

        msg('finished!', level=2, silent=silent)

        msg('first {0} eigenvalues:'.format(self.num_eigvalues_print), level=1,
            silent=silent)
        for eig in eigvals[:self.num_eigvalues_print]:
            msg('{0}'.format(eig), level=2, silent=silent)
        self.analysis.last_analysis = 'lb'


    def freq(self, atype=4, tol=0, sparse_solver=True, silent=False,
             sort=True, damping=False, reduced_dof=False):
        """Performs a natural frequency analysis

        .. note:: This will be deprecated soon, use
                  :func:`.compmech.analysis.freq`.

        The following parameters of the will affect the linear buckling
        analysis:

        =======================    =====================================
        parameter                  description
        =======================    =====================================
        ``num_eigvalues``        Number of eigenvalues to be extracted
        ``num_eigvalues_print``    Number of eigenvalues to print after
                                   the analysis is completed
        =======================    =====================================

        Parameters
        ----------
        atype : int, optional
            Tells which analysis type should be performed:

            - ``1`` : considers k0, kA and kG0 (and cA depending on 'damping')
            - ``2`` : considers k0 and kA (and cA depending on 'damping')
            - ``3`` : considers k0 and kG0
            - ``4`` : considers k0 only

        tol : float, optional
            A tolerance value passed to ``scipy.sparse.linalg.eigs``.
        sparse_solver : bool, optional
            Tells if solver :func:`scipy.linalg.eig` or
            :func:`scipy.sparse.linalg.eigs` should be used.

            .. note:: It is recommended ``sparse_solver=False``, because it
                      was verified that the sparse solver becomes unstable
                      for some cases, though the sparse solver is faster.
        silent : bool, optional
            A boolean to tell whether the log messages should be printed.
        sort : bool, optional
            Sort the output eigenvalues and eigenmodes.
        damping : bool, optinal
            If aerodynamic damping should be taken into account.
        reduced_dof : bool, optional
            Considers only the contributions of `v` and `w` to the stiffness
            matrix and accelerates the run. Only effective when
            ``sparse_solver=False``.

        Notes
        -----
        The extracted eigenvalues are stored in the ``eigvals`` parameter and
        the `i^{th}` eigenvector in the ``eigvecs[:, i-1]`` parameter.

        """
        msg('Running frequency analysis...', silent=silent)

        self.calc_k0(silent=silent)
        self.calc_kM(silent=silent)
        if atype == 1:
            self.calc_kG0(silent=silent)
            self.calc_kA(silent=silent)
            if damping:
                self.calc_cA(silent=silent)
                K = self.k0 + self.kA + self.kG0 + self.cA
            else:
                K = self.k0 + self.kA + self.kG0

        elif atype == 2:
            self.calc_kA(silent=silent)
            K = self.k0 + self.kA
            if damping:
                self.calc_cA(silent=silent)
                K = self.k0 + self.kA + self.cA
            else:
                K = self.k0 + self.kA

        elif atype == 3:
            self.calc_kG0(silent=silent)
            K = self.k0 + self.kG0

        elif atype == 4:
            K = self.k0

        M = self.kM

        msg('Eigenvalue solver... ', level=2, silent=silent)
        k = min(self.num_eigvalues, M.shape[0]-2)
        if sparse_solver:
            if damping:
                raise NotImplementedError('Damping with sparse_solver not implemented!')

            msg('eigs() solver...', level=3, silent=silent)
            sizebkp = M.shape[0]
            K, M, used_cols = remove_null_cols(K, M, silent=silent,
                    level=3)
            eigvals, peigvecs = eigs(A=K, k=k, which='LM', M=M, tol=tol,
                                     sigma=-1.)
            eigvecs = np.zeros((sizebkp, self.num_eigvalues),
                               dtype=peigvecs.dtype)
            eigvecs[used_cols, :] = peigvecs

            eigvals = np.sqrt(eigvals) # omega^2 to omega, in rad/s

        else:
            msg('eig() solver...', level=3, silent=silent)
            M = M.toarray()
            K = K.toarray()
            sizebkp = M.shape[0]
            col_sum = M.sum(axis=0)
            check = col_sum != 0
            used_cols = np.arange(M.shape[0])[check]
            M = M[:, check][check, :]
            K = K[:, check][check, :]

            if reduced_dof:
                i = np.arange(M.shape[0])
                take = np.column_stack((i[1::3], i[2::3])).flatten()
                M = M[:, take][take, :]
                K = K[:, take][take, :]
            if not damping:
                M = -M
            else:
                size = M.shape[0]
                cA = self.cA.toarray()
                cA = cA[:, check][check, :]
                if reduced_dof:
                    cA = cA[:, take][take, :]
                I = np.identity(M.shape[0])
                Z = np.zeros_like(M)
                M = np.row_stack((np.column_stack((I, Z)),
                                  np.column_stack((Z, -M))))
                K = np.row_stack((np.column_stack((Z, -I)),
                                  np.column_stack((K, cA))))

            eigvals, peigvecs = eig(a=M, b=K)

            eigvecs = np.zeros((sizebkp, K.shape[0]),
                               dtype=peigvecs.dtype)
            eigvecs[check, :] = peigvecs

            if not damping:
                eigvals = np.sqrt(-1./eigvals) # -1/omega^2 to omega, in rad/s
                eigvals = eigvals
            else:
                eigvals = -1./eigvals # -1/omega to omega, in rad/s
                shape = eigvals.shape
                eigvals = eigvals[:shape[0]//2]
                eigvecs = eigvecs[:eigvecs.shape[0]//2, :shape[0]//2]

        msg('finished!', level=3, silent=silent)

        if sort:
            if damping:
                higher_zero = eigvals.real > 1e-6

                eigvals = eigvals[higher_zero]
                eigvecs = eigvecs[:, higher_zero]

                sort_ind = np.lexsort((np.round(eigvals.imag, 1),
                                       np.round(eigvals.real, 0)))
                eigvals = eigvals[sort_ind]
                eigvecs = eigvecs[:, sort_ind]

            else:
                sort_ind = np.lexsort((np.round(eigvals.imag, 1),
                                       np.round(eigvals.real, 1)))
                eigvals = eigvals[sort_ind]
                eigvecs = eigvecs[:, sort_ind]

                higher_zero = eigvals.real > 1e-6

                eigvals = eigvals[higher_zero]
                eigvecs = eigvecs[:, higher_zero]

        if not sparse_solver and reduced_dof:
            new_eigvecs = np.zeros((3*eigvecs.shape[0]//2, eigvecs.shape[1]),
                    dtype=eigvecs.dtype)
            new_eigvecs[take, :] = eigvecs
            eigvecs = new_eigvecs

        self.eigvals = eigvals
        self.eigvecs = eigvecs

        msg('finished!', level=2, silent=silent)

        msg('first {0} eigenvalues:'.format(self.num_eigvalues_print), level=1,
            silent=silent)
        for eigval in eigvals[:self.num_eigvalues_print]:
            msg('{0} rad/s'.format(eigval), level=2, silent=silent)
        self.analysis.last_analysis = 'freq'


    def uvw(self, c, xs=None, ys=None, gridx=300, gridy=300):
        r"""Calculates the displacement field

        For a given full set of Ritz constants ``c``, the displacement
        field is calculated and stored in the parameters
        ``u``, ``v``, ``w``, ``phix``, ``phiy`` of the ``Panel`` object.

        Parameters
        ----------
        c : float
            The full set of Ritz constants
        xs : np.ndarray
            The `x` positions where to calculate the displacement field.
            Default is ``None`` and the method ``_default_field`` is used.
        ys : np.ndarray
            The ``y`` positions where to calculate the displacement field.
            Default is ``None`` and the method ``_default_field`` is used.
        gridx : int
            Number of points along the `x` axis where to calculate the
            displacement field.
        gridy : int
            Number of points along the `y` where to calculate the
            displacement field.

        Returns
        -------
        out : tuple
            A tuple of ``np.ndarrays`` containing
            ``(u, v, w, phix, phiy)``.

        Notes
        -----
        The returned values ``u```, ``v``, ``w``, ``phix``, ``phiy`` are
        stored as parameters with the same name in the ``Panel`` object.

        """
        c = np.ascontiguousarray(c, dtype=DOUBLE)

        xs, ys, xshape, tshape = self._default_field(xs, ys, gridx, gridy)
        fuvw = modelDB.db[self.model]['field'].fuvw
        us, vs, ws, phixs, phiys = fuvw(c, self, xs, ys, self.out_num_cores)

        self.u = us.reshape(xshape)
        self.v = vs.reshape(xshape)
        self.w = ws.reshape(xshape)
        self.phix = phixs.reshape(xshape)
        self.phiy = phiys.reshape(xshape)

        return self.u, self.v, self.w, self.phix, self.phiy


    def strain(self, c, xs=None, ys=None, gridx=300, gridy=300):
        r"""Calculates the strain field

        Parameters
        ----------
        c : np.ndarray
            The Ritz constants vector to be used for the strain field
            calculation.
        xs : np.ndarray, optional
            The `x` coordinates where to calculate the strains.
        ys : np.ndarray, optional
            The `y` coordinates where to calculate the strains, must
            have the same shape as ``xs``.
        gridx : int, optional
            When ``xs`` and ``ys`` are not supplied, ``gridx`` and ``gridy``
            are used.
        gridy : int, optional
            When ``xs`` and ``ys`` are not supplied, ``gridx`` and ``gridy``
            are used.

        """
        c = np.ascontiguousarray(c, dtype=DOUBLE)

        xs, ys, xshape, tshape = self._default_field(xs, ys, gridx, gridy)

        model = self.model
        NL_kinematics = model.split('_')[1]
        fstrain = modelDB.db[model]['field'].fstrain
        e_num = modelDB.db[model]['e_num']

        if 'donnell' in NL_kinematics:
            int_NL_kinematics = 0
        elif 'sanders' in NL_kinematics:
            int_NL_kinematics = 1
        else:
            raise NotImplementedError(
                '{} is not a valid NL_kinematics option'.format(NL_kinematics))

        es = fstrain(c, xs, ys, self.a, self.b, self.m, self.n, self.c0,
                self.m0, self.n0, self.funcnum, int_NL_kinematics,
                self.out_num_cores)

        return es.reshape((xshape + (e_num,)))


    def stress(self, c, xs=None, ys=None, gridx=300, gridy=300):
        r"""Calculates the stress field

        Parameters
        ----------
        c : np.ndarray
            The Ritz constants vector to be used for the strain field
            calculation.
        xs : np.ndarray, optional
            The `x` coordinates where to calculate the strains.
        ys : np.ndarray, optional
            The `y` coordinates where to calculate the strains, must
            have the same shape as ``xs``.
        gridx : int, optional
            When ``xs`` and ``ys`` are not supplied, ``gridx`` and ``gridy``
            are used.
        gridy : int, optional
            When ``xs`` and ``ys`` are not supplied, ``gridx`` and ``gridy``
            are used.

        """
        c = np.ascontiguousarray(c, dtype=DOUBLE)

        xs, ys, xshape, tshape = self._default_field(xs, ys, gridx, gridy)

        model = self.model
        NL_kinematics = model.split('_')[1]
        fstress = modelDB.db[model]['field'].fstress
        e_num = modelDB.db[model]['e_num']

        if 'donnell' in NL_kinematics:
            int_NL_kinematics = 0
        elif 'sanders' in NL_kinematics:
            int_NL_kinematics = 1
        else:
            raise NotImplementedError(
                    '{} is not a valid NL_kinematics option'.format(
                    NL_kinematics))
        Ns = fstress(c, self.F, xs, ys, self.a, self.b, self.m, self.n,
                self.c0, self.m0, self.n0,
                self.funcnum, int_NL_kinematics, self.out_num_cores)

        return Ns.reshape((xshape + (e_num,)))


    def add_force(self, x, y, fx, fy, fz, cte=True):
        r"""Add a punctual force with three components

        Parameters
        ----------
        x : float
            The `x` position.
        y : float
            The `y` position in radians.
        fx : float
            The `x` component of the force vector.
        fy : float
            The `y` component of the force vector.
        fz : float
            The `z` component of the force vector.
        cte : bool, optional
            Constant forces are not incremented during the non-linear
            analysis.

        """
        if cte:
            self.forces.append([x, y, fx, fy, fz])
        else:
            self.forces_inc.append([x, y, fx, fy, fz])


    def calc_fext(self, inc=1., silent=False):
        """Calculates the external force vector `\{F_{ext}\}`

        Recall that:

        .. math::

            \{F_{ext}\}=\{{F_{ext}}_0\} + \{{F_{ext}}_\lambda\}

        such that the terms in `\{{F_{ext}}_0\}` are constant and the terms in
        `\{{F_{ext}}_\lambda\}` will be scaled by the parameter ``inc``.

        Parameters
        ----------
        inc : float, optional
            Since this function is called during the non-linear analysis,
            ``inc`` will multiply the terms `\{{F_{ext}}_\lambda\}`.

        silent : bool, optional
            A boolean to tell whether the log messages should be printed.

        Returns
        -------
        fext : np.ndarray
            The external force vector

        """
        self._rebuild()
        msg('Calculating external forces...', level=2, silent=silent)

        model = self.model
        if not model in modelDB.db.keys():
            raise ValueError(
                    '{} is not a valid model option'.format(model))
        db = modelDB.db
        num = db[model]['num']
        dofs = db[model]['dofs']
        fg = db[model]['field'].fg

        size = self.get_size()
        g = np.zeros((dofs, size), dtype=DOUBLE)
        fext = np.zeros(size, dtype=DOUBLE)

        # non-incrementable punctual forces
        for i, force in enumerate(self.forces):
            x, y, fx, fy, fz = force
            fg(g, x, y, self)
            if dofs == 3:
                fpt = np.array([[fx, fy, fz]])
            elif dofs == 5:
                fpt = np.array([[fx, fy, fz, 0, 0]])
            fext += fpt.dot(g).ravel()

        # incrementable punctual forces
        for i, force in enumerate(self.forces_inc):
            x, y, fx, fy, fz = force
            fg(g, x, y, self)
            if dofs == 3: #CLT
                fpt = np.array([[fx, fy, fz]])*inc
            elif dofs == 5: #FSDT
                fpt = np.array([[fx, fy, fz, 0, 0]])*inc
            fext += fpt.dot(g).ravel()

        return fext


    def calc_fint(self, c, size=None, col0=0, silent=False, nx=None,
            ny=None, Fnxny=None, inc=None):
        #TODO inc not needed here; otherwise with prescribed displacements
        msg('Calculating internal forces...', level=2, silent=silent)
        model = self.model
        if not model in modelDB.db.keys():
            raise ValueError(
                    '{0} is not a valid model option'.format(model))
        matrices_num = modelDB.db[model].get('matrices_num')
        if matrices_num is None:
            raise ValuError('matrices_num not implemented for model {0}'.
                    format(model))
        calc_fint = getattr(matrices_num, 'calc_fint', None)
        if calc_fint is None:
            raise ValuError('calc_fint not implemented for model {0}'.
                    format(model))

        if size is None:
            size = self.get_size()

        alphadeg = self.alphadeg if self.alphadeg is not None else 0.
        self.alpharad = deg2rad(alphadeg)
        self.r = self.r if self.r is not None else 0.
        nx = self.nx if nx is None else nx
        ny = self.ny if ny is None else ny
        Fnxny = self.F if Fnxny is None else Fnxny

        c = np.ascontiguousarray(c, dtype=DOUBLE)
        fint = calc_fint(c, Fnxny, self, size, col0, nx, ny)

        gc.collect()

        msg('finished!', level=2, silent=silent)

        return fint


    def static(self, NLgeom=False, silent=False):
        """Static analysis for cones and cylinders

        .. note:: This will be deprecated soon, use
                  :func:`.compmech.analysis.static`.

        The analysis can be linear or geometrically non-linear. See
        :class:`.Analysis` for further details about the parameters
        controlling the non-linear analysis.

        Parameters
        ----------
        NLgeom : bool
            Flag to indicate whether a linear or a non-linear analysis is to
            be performed.
        silent : bool, optional
            A boolean to tell whether the log messages should be printed.

        Returns
        -------
        cs : list
            A list containing the Ritz constants for each load increment of
            the static analysis. The list will have only one entry in case
            of a linear analysis.

        Notes
        -----
        The returned ``cs`` is stored in ``self.analysis.cs``. The actual
        increments used in the non-linear analysis are stored in the
        ``self.analysis.increments`` parameter.

        """
        self._rebuild()
        self.analysis.line_search = False
        self.analysis.kT_initial_state = False
        self.analysis.compute_every_n = 6

        if NLgeom and not modelDB.db[self.model]['non-linear static']:
            msg('________________________________________________',
                silent=silent)
            msg('', silent=silent)
            warn('Model {} cannot be used in non-linear static analysis!'.
                 format(self.model), silent=silent)
            msg('________________________________________________',
                silent=silent)
            raise
        elif not NLgeom and not modelDB.db[self.model]['matrices']:
            msg('________________________________________________',
                level=1, silent=silent)
            msg('', level=1, silent=silent)
            warn('Model {} cannot be used in linear static analysis!'.
                 format(self.model), level=1, silent=silent)
            msg('________________________________________________',
                level=1, silent=silent)
            raise
        self.analysis.static(NLgeom=NLgeom, silent=silent)
        self.increments = self.analysis.increments

        return self.analysis.cs


    def plot(self, c, invert_y=False, vec='w',
             deform_u=False, deform_u_sf=100.,
             filename='',
             ax=None, figsize=(3.5, 2.), save=True,
             title='',
             colorbar=False, cbar_nticks=2, cbar_format=None,
             cbar_title='', cbar_fontsize=10,
             aspect='equal', clean=True, dpi=400,
             texts=[], xs=None, ys=None, gridx=300, gridy=300,
             num_levels=400, vecmin=None, vecmax=None, plotoffsetxs=0.,
             plotoffsetys=0.):
        r"""Contour plot for a Ritz constants vector.

        Parameters
        ----------
        c : np.ndarray
            The Ritz constants that will be used to compute the field contour.
        vec : str, optional
            Can be one of the components:

            - Displacement: ``'u'``, ``'v'``, ``'w'``, ``'phix'``, ``'phiy'``
            - Strain: ``'exx'``, ``'eyy'``, ``'gxy'``, ``'kxx'``, ``'kyy'``,
              ``'kxy'``, ``'gyz'``, ``'gxz'``
            - Stress: ``'Nxx'``, ``'Nyy'``, ``'Nxy'``, ``'Mxx'``, ``'Myy'``,
              ``'Mxy'``, ``'Qy'``, ``'Qx'``
        deform_u : bool, optional
            If ``True`` the contour plot will look deformed.
        deform_u_sf : float, optional
            The scaling factor used to deform the contour.
        invert_y : bool, optional
            Inverts the `y` axis of the plot.
        save : bool, optional
            Flag telling whether the contour should be saved to an image file.
        dpi : int, optional
            Resolution of the saved file in dots per inch.
        filename : str, optional
            The file name for the generated image file. If no value is given,
            the `name` parameter of the ``Panel`` object will be used.
        ax : AxesSubplot, optional
            When ``ax`` is given, the contour plot will be created inside it.
        figsize : tuple, optional
            The figure size given by ``(width, height)``.
        title : str, optional
            If any string is given it is added as title to the contour plot.
        colorbar : bool, optional
            If a colorbar should be added to the contour plot.
        cbar_nticks : int, optional
            Number of ticks added to the colorbar.
        cbar_format : [ None | format string | Formatter object ], optional
            See the ``matplotlib.pyplot.colorbar`` documentation.
        cbar_fontsize : int, optional
            Fontsize of the colorbar labels.
        cbar_title : str, optional
            Colorbar title. If ``cbar_title == ''`` no title is added.
        aspect : str, optional
            String that will be passed to the ``AxesSubplot.set_aspect()``
            method.
        clean : bool, optional
            Clean axes ticks, grids, spines etc.
        xs : np.ndarray, optional
            The `x` positions where to calculate the displacement field.
            Default is ``None`` and the method ``_default_field`` is used.
        ys : np.ndarray, optional
            The ``y`` positions where to calculate the displacement field.
            Default is ``None`` and the method ``_default_field`` is used.
        gridx : int, optional
            Number of points along the `x` axis where to calculate the
            displacement field.
        gridy : int, optional
            Number of points along the `y` where to calculate the
            displacement field.
        num_levels : int, optional
            Number of contour levels (higher values make the contour smoother).
        vecmin : float, optional
            Minimum value for the contour scale (useful to compare with other
            results). If not specified it will be taken from the calculated
            field.
        vecmax : float, optional
            Maximum value for the contour scale.

        Returns
        -------
        ax : matplotlib.axes.Axes
            The Matplotlib object that can be used to modify the current plot
            if needed.

        """
        msg('Plotting contour...')

        ubkp, vbkp, wbkp, phixbkp, phiybkp = (self.u, self.v, self.w,
                                              self.phix, self.phiy)

        import matplotlib
        if platform.system().lower() == 'linux':
            matplotlib.use('Agg')
        import matplotlib.pyplot as plt

        msg('Computing field variables...', level=1)
        displs = ['u', 'v', 'w', 'phix', 'phiy']
        strains = ['exx', 'eyy', 'gxy', 'kxx', 'kyy', 'kxy', 'gyz', 'gxz']
        stresses = ['Nxx', 'Nyy', 'Nxy', 'Mxx', 'Myy', 'Mxy', 'Qy', 'Qx']
        if vec in displs:
            self.uvw(c, xs=xs, ys=ys, gridx=gridx, gridy=gridy)
            field = getattr(self, vec)
        elif vec in strains:
            es = self.strain(c, xs=xs, ys=ys,
                             gridx=gridx, gridy=gridy)
            field = es[..., strains.index(vec)]
        elif vec in stresses:
            Ns = self.stress(c, xs=xs, ys=ys,
                             gridx=gridx, gridy=gridy)
            field = Ns[..., stresses.index(vec)]
        else:
            raise ValueError(
                    '{0} is not a valid vec parameter value!'.format(vec))
        msg('Finished!', level=1)

        Xs = self.Xs + plotoffsetxs
        Ys = self.Ys + plotoffsetys

        if vecmin is None:
            vecmin = field.min()
        if vecmax is None:
            vecmax = field.max()

        levels = linspace(vecmin, vecmax, num_levels)

        if ax is None:
            fig = plt.figure(figsize=figsize)
            ax = fig.add_subplot(111)
        else:
            if isinstance(ax, matplotlib.axes.Axes):
                ax = ax
                fig = ax.figure
                save = False
            else:
                raise ValueError('ax must be an Axes object')

        x = Ys # in matplotlib x goes vertically (axis=0)
        y = Xs # and y goes horizontally (axis=1)

        if deform_u:
            if vec in displs:
                pass
            else:
                self.uvw(c, xs=xs, ys=ys, gridx=gridx, gridy=gridy)
            field_u = self.u
            field_v = self.v
            y -= deform_u_sf*field_u
            x += deform_u_sf*field_v
        contour = ax.contourf(x, y, field, levels=levels)

        if colorbar:
            from mpl_toolkits.axes_grid1 import make_axes_locatable

            fsize = cbar_fontsize
            divider = make_axes_locatable(ax)
            cax = divider.append_axes('right', size='5%', pad=0.05)
            cbarticks = linspace(vecmin, vecmax, cbar_nticks)
            cbar = plt.colorbar(contour, ticks=cbarticks, format=cbar_format,
                                cax=cax)
            if cbar_title:
                cax.text(0.5, 1.05, cbar_title, horizontalalignment='center',
                         verticalalignment='bottom', fontsize=fsize)
            cbar.outline.remove()
            cbar.ax.tick_params(labelsize=fsize, pad=0., tick2On=False)

        if invert_y == True:
            ax.invert_yaxis()
        ax.invert_xaxis()

        if title != '':
            ax.set_title(str(title))

        fig.tight_layout()
        ax.set_aspect(aspect)

        ax.grid(False)
        ax.set_frame_on(False)
        if clean:
            ax.xaxis.set_ticks_position('none')
            ax.yaxis.set_ticks_position('none')
            ax.xaxis.set_ticklabels([])
            ax.yaxis.set_ticklabels([])
        else:
            ax.xaxis.set_ticks_position('bottom')
            ax.yaxis.set_ticks_position('left')

        for kwargs in texts:
            ax.text(transform=ax.transAxes, **kwargs)

        if save:
            if not filename:
                filename = 'test.png'
            fig.savefig(filename, transparent=True,
                        bbox_inches='tight', pad_inches=0.05, dpi=dpi)
            plt.close()

        if ubkp is not None:
            self.u = ubkp
        if vbkp is not None:
            self.v = vbkp
        if wbkp is not None:
            self.w = wbkp
        if phixbkp is not None:
            self.phix = phixbkp
        if phiybkp is not None:
            self.phiy = phiybkp

        msg('finished!')

        return ax


    def save(self):
        """Save the ``Panel`` object using ``pickle``

        Notes
        -----
        The pickled file will have the name stored in ``Panel.name``
        followed by a ``'.Panel'`` extension.

        """
        name = self.name + '.Panel'
        msg('Saving Panel to {}'.format(name))

        self._clear_matrices()

        with open(name, 'wb') as f:
            pickle.dump(self, f, protocol=pickle.HIGHEST_PROTOCOL)

