from __future__ import division
import gc
import os
import sys
import traceback
from collections import Iterable
import time
import cPickle
import __main__

import numpy as np
from scipy.sparse import coo_matrix, csr_matrix, csc_matrix
from scipy.sparse.linalg import eigsh
from scipy.optimize import leastsq
from numpy import linspace, pi, cos, sin, tan, deg2rad

from conecylDB import ccs, laminaprops
import compmech.composite.laminate as laminate
from compmech.analysis import Analysis
from compmech.logger import msg, warn, error
from compmech.sparse import remove_null_cols, make_symmetric
from compmech.constants import DOUBLE
import modelDB
from modelDB import get_model


def load(name):
    if '.ConeCyl' in name:
        cc = cPickle.load(open(name, 'rb'))
    else:
        cc = cPickle.load(open(name + '.ConeCyl', 'rb'))
    cc.analysis.calc_fext = cc.calc_fext
    cc.analysis.calc_k0 = cc.calc_k0
    cc.analysis.calc_fint = cc.calc_fint
    cc.analysis.calc_kT = cc.calc_kT
    return cc


class ConeCyl(object):
    """
    """
    __slots__ = ['_load_rebuilt', 'name', 'alphadeg', 'alpharad', 'r1', 'r2',
            'L', 'H', 'h', 'K', 'is_cylinder', 'inf', 'zero',
            'bc', 'kuBot', 'kvBot', 'kwBot', 'kphixBot', 'kphitBot', 'kuTop',
            'kvTop', 'kwTop', 'kphixTop', 'kphitTop', 'model', 'm1', 'm2',
            'size', 'n2', 's', 'nx', 'nt', 'ni_num_cores', 'ni_method',
            'forces', 'forces_inc', 'P',
            'P_inc', 'pdC', 'Fc', 'Nxxtop', 'uTM', 'c0', 'm0', 'n0',
            'funcnum', 'pdT', 'T', 'T_inc', 'thetaTdeg', 'thetaTrad', 'pdLA',
            'tLAdeg', 'tLArad', 'betadeg', 'betarad', 'xiLA', 'MLA', 'LA',
            'num0', 'excluded_dofs', 'excluded_dofs_ck', 'sina', 'cosa',
            'laminaprop', 'plyt', 'laminaprops', 'stack', 'plyts', 'F',
            'F_reuse', 'force_orthotropic_laminate', 'E11', 'nu',
            'num_eigvalues', 'num_eigvalues_print',

            'analysis', 'with_k0L', 'with_kLL',
            'cs', 'increments', 'outputs',
            'eigvals', 'eigvecs',

            'k0', 'k0uk', 'k0uu', 'kTuk', 'kTuu', 'kG0', 'kG0_Fc', 'kG0_P',
            'kG0_T', 'kG', 'kGuu', 'kL', 'kLuu', 'lam', 'u', 'v', 'w', 'phix',
            'phit', 'Xs', 'Ts',

            'out_num_cores',
        ]

    def __init__(self):
        self.name = ''

        # geometry
        self.alphadeg = 0.
        self.alpharad = 0.
        self.r1 = None
        self.r2 = None
        self.L = None
        self.H = None
        self.h = None # total thickness, required for isotropic shells
        self.is_cylinder = None

        # boundary conditions
        self.inf = 1.e8 # used to define high stiffnesses
        self.zero = 0. # used to define zero stiffnesses
        self.bc = None
        self.kuBot = self.inf
        self.kvBot = self.inf
        self.kwBot = self.inf
        self.kphixBot = 0.
        self.kphitBot = 0.
        self.kuTop = self.inf
        self.kvTop = self.inf
        self.kwTop = self.inf
        self.kphixTop = 0.
        self.kphitTop = 0.

        # default equations
        self.model = 'clpt_donnell_bc1'

        # approximation series
        self.m1 = 120
        self.m2 = 25
        self.n2 = 45

        # analytical integration for cones
        self.s = 79

        # numerical integration
        self.nx = 120
        self.nt = 180
        self.ni_num_cores = 4
        self.ni_method = 'trapz2d'

        # punctual loads
        self.forces = []
        self.forces_inc = []

        # internal pressure measured in force/area
        self.P = 0.
        self.P_inc = 0.

        # axial compression
        self.pdC = False
        self.Fc = None
        self.Nxxtop = None
        self.uTM = 0.
        self._load_rebuilt = False

        # initial imperfection
        self.c0 = None
        self.m0 = 0
        self.n0 = 0
        self.funcnum = 2

        # torsion
        self.pdT = True
        self.T = 0.
        self.T_inc = 0.
        self.thetaTdeg = 0.
        self.thetaTrad = 0.

        # load asymmetry (la)
        self.pdLA = True
        self.tLAdeg = 0.
        self.tLArad = 0.
        self.betadeg = 0.
        self.betarad = 0.
        self.xiLA = None
        self.MLA = None
        self.LA = None

        self.num0 = 3
        self.excluded_dofs = []
        self.excluded_dofs_ck = []

        self.sina = None
        self.cosa = None

        # material
        self.laminaprop = None
        self.plyt = None
        self.laminaprops = []
        self.stack = []
        self.plyts = []

        # constitutive law
        self.F_reuse = None
        self.F = None
        self.force_orthotropic_laminate = False
        self.E11 = None
        self.nu = None
        self.K = 5/6.

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

        # output queries
        self.out_num_cores = 4
        self.cs = []
        self.increments = []

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

        # outputs
        self.outputs = {}

        self._clear_matrices()


    def _clear_matrices(self):
        self.k0 = None
        self.k0uk = None
        self.k0uu = None
        self.kTuk = None
        self.kTuu = None
        self.kG0 = None
        self.kG0_Fc = None
        self.kG0_P = None
        self.kG0_T = None
        self.kG = None
        self.kGuu = None
        self.kL = None
        self.kLuu = None
        self.lam = None
        self.u = None
        self.v = None
        self.w = None
        self.phix = None
        self.phit = None
        self.Xs = None
        self.Ts = None
        self.Nxxtop = None

        gc.collect()


    def _rebuild(self):
        if not self.name:
            try:
                self.name = os.path.basename(__main__.__file__).split('.py')[0]
            except AttributeError:
                warn('ConeCyl name unchanged')

        if self.k0 is not None:
            if self.k0.shape[0] != self.get_size():
                self._clear_matrices()
                self._load_rebuilt = False
                self._rebuild()

        self.model = self.model.lower()
        model_dict = get_model(self.model)

        # boundary conditions
        inf = self.inf
        zero = self.zero

        if inf > 1.e8:
            warn('"inf" parameter reduced to 1.e8 due to the verified ' +
                 'numerical instability for higher values', level=2)
            inf = 1.e8

        if self.bc is not None:
            bc = self.bc.lower()

            if '_' in bc:
                # different bc for Bot and Top
                bc_Bot, bc_Top = self.bc.split('_')
            elif '-' in bc:
                # different bc for Bot and Top
                bc_Bot, bc_Top = self.bc.split('-')
            else:
                bc_Bot = bc_Top = bc

            bcs = dict(bc_Bot=bc_Bot, bc_Top=bc_Top)
            for k in bcs.keys():
                sufix = k.split('_')[1] # Bot or Top
                if bcs[k] == 'ss1':
                    setattr(self, 'ku' + sufix, inf)
                    setattr(self, 'kv' + sufix, inf)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, zero)
                    setattr(self, 'kphit' + sufix, zero)
                elif bcs[k] == 'ss2':
                    setattr(self, 'ku' + sufix, zero)
                    setattr(self, 'kv' + sufix, inf)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, zero)
                    setattr(self, 'kphit' + sufix, zero)
                elif bcs[k] == 'ss3':
                    setattr(self, 'ku' + sufix, inf)
                    setattr(self, 'kv' + sufix, zero)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, zero)
                    setattr(self, 'kphit' + sufix, zero)
                elif bcs[k] == 'ss4':
                    setattr(self, 'ku' + sufix, zero)
                    setattr(self, 'kv' + sufix, zero)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, zero)
                    setattr(self, 'kphit' + sufix, zero)

                elif bcs[k] == 'cc1':
                    setattr(self, 'ku' + sufix, inf)
                    setattr(self, 'kv' + sufix, inf)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, inf)
                    setattr(self, 'kphit' + sufix, zero)
                elif bcs[k] == 'cc2':
                    setattr(self, 'ku' + sufix, zero)
                    setattr(self, 'kv' + sufix, inf)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, inf)
                    setattr(self, 'kphit' + sufix, zero)
                elif bcs[k] == 'cc3':
                    setattr(self, 'ku' + sufix, inf)
                    setattr(self, 'kv' + sufix, zero)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, inf)
                    setattr(self, 'kphit' + sufix, zero)
                elif bcs[k] == 'cc4':
                    setattr(self, 'ku' + sufix, zero)
                    setattr(self, 'kv' + sufix, zero)
                    setattr(self, 'kw' + sufix, inf)
                    setattr(self, 'kphix' + sufix, inf)
                    setattr(self, 'kphit' + sufix, zero)

                elif bcs[k] == 'free':
                    setattr(self, 'ku' + sufix, zero)
                    setattr(self, 'kv' + sufix, zero)
                    setattr(self, 'kw' + sufix, zero)
                    setattr(self, 'kphix' + sufix, zero)
                    setattr(self, 'kphit' + sufix, zero)
                else:
                    text = '"{}" is not a valid boundary condition!'.format(bc)
                    raise ValueError(text)

        self.alpharad = deg2rad(self.alphadeg)
        self.sina = sin(self.alpharad)
        self.cosa = cos(self.alpharad)

        if not self.H and not self.L:
            self.H = (self.r1-self.r2)/tan(self.alpharad)
        if self.H and not self.L:
            self.L = self.H/self.cosa
        if self.L and not self.H:
            self.H = self.L*self.cosa

        if not self.r2:
            if not self.r1:
                raise ValueError('Radius "r1" or "r2" must be specified')
            else:
                self.r2 = self.r1 - self.L*self.sina
        else:
            self.r1 = self.r2 + self.L*self.sina

        self.thetaTrad = deg2rad(self.thetaTdeg)

        self.tLArad = deg2rad(self.tLAdeg)

        self.betarad = deg2rad(self.betadeg)
        self.LA = self.r2*tan(self.betarad)

        if not self.laminaprops:
            self.laminaprops = [self.laminaprop for i in self.stack]
        if not self.plyts:
            self.plyts = [self.plyt for i in self.stack]

        if self.alpharad == 0:
            self.is_cylinder = True
        else:
            self.is_cylinder = False

        self.excluded_dofs = []
        self.excluded_dofs_ck = []
        if self.pdC:
            self.excluded_dofs.append(0)
            self.excluded_dofs_ck.append(self.uTM)
        if self.pdT:
            self.excluded_dofs.append(1)
            self.excluded_dofs_ck.append(self.thetaTrad)
        if self.pdLA:
            self.excluded_dofs.append(2)
            self.excluded_dofs_ck.append(self.LA)
        else:
            raise NotImplementedError('pdLA == False is giving wrong results!')

        if self.nx < 4*self.m2:
            warn('Number of integration points along x too small')
        if self.nt < 4*self.n2:
            warn('Number of integration points along theta too small')

        if self.laminaprop is None:
            h = self.h
            E11 = self.E11
            nu = self.nu
            if h is None or E11 is None or nu is None:
                raise ValueError(
                        'laminaprop or (E11, nu and h) must be defined')

            G12 = E11/(2*(1 + nu))
            A11 = E11*h/(1 - nu**2)
            A12 = nu*E11*h/(1 - nu**2)
            A16 = 0
            A22 = E11*h/(1 - nu**2)
            A26 = 0
            A66 = G12*h
            D11 = E11*h**3/(12*(1 - nu**2))
            D12 = nu*E11*h**3/(12*(1 - nu**2))
            D16 = 0
            D22 = E11*h**3/(12*(1 - nu**2))
            D26 = 0
            D66 = G12*h**3/12

            # TODO, what if FSDT is used?
            if 'fsdt' in self.model:
                raise NotImplementedError(
                        'For FSDT laminaprop must be defined!')

            self.F = np.array([[A11, A12, A16, 0, 0, 0],
                               [A12, A22, A26, 0, 0, 0],
                               [A16, A26, A66, 0, 0, 0],
                               [0, 0, 0, D11, D12, D16],
                               [0, 0, 0, D12, D22, D26],
                               [0, 0, 0, D16, D26, D66]])

        if self.c0 is not None:
            self.analysis.kT_initial_state = True

        if self.Nxxtop is not None and self._load_rebuilt:
            return

        # axial load
        if self.Nxxtop is not None:
            if type(self.Nxxtop) in (int, float):
                Nxxtop0 = self.Nxxtop
                self.Nxxtop = np.zeros(2*self.n2+1, dtype=DOUBLE)
                self.Nxxtop[0] = Nxxtop0

            check = False
            if isinstance(self.Nxxtop, np.ndarray):
                if self.Nxxtop.ndim == 1:
                    assert self.Nxxtop.shape[0] == (2*self.n2+1)
                    check = True
            if not check:
                raise ValueError('Invalid Nxxtop input')

        else:
            self.Nxxtop = np.zeros(2*self.n2+1, dtype=DOUBLE)

        if self.Fc is not None:
            self.Nxxtop[0] = self.Fc/(2*pi*self.r2*self.cosa)
            msg('Nxxtop[0] calculated from Fc', level=2)
            if self.MLA is None:
                if self.xiLA is not None:
                    self.MLA = self.xiLA*self.Fc
                    msg('MLA calculated from xiLA', level=2)

        if self.MLA is not None:
            self.Nxxtop[2] = self.MLA/(pi*self.r2**2*self.cosa)
            msg('Nxxtop[2] calculated from MLA', level=2)

        self._load_rebuilt = True


    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.

        """
        model_dict = get_model(self.model)
        num0 = model_dict['num0']
        num1 = model_dict['num1']
        num2 = model_dict['num2']
        self.size = num0 + num1*self.m1 + num2*self.m2*self.n2
        return self.size


    def from_DB(self, name):
        """Load cone / cylinder data from the local database

        Parameters
        ----------
        name : str
            A key contained in the ``ccs`` dictionary of module
            :mod:`compmech.conecyl.conecylDB`.

        """
        try:
            attrs = ['r1', 'r2', 'H', 'L', 'alphadeg', 'plyt', 'stack']
            cc = ccs[name]
            self.laminaprop = laminaprops[cc['laminapropKey']]
            for attr in attrs:
                setattr(self, attr, cc.get(attr, getattr(self, attr)))
        except:
            raise ValueError('Invalid data-base entry!')


    def exclude_dofs_matrix(self, k, return_kkk=False,
                                     return_kku=False,
                                     return_kuk=False):
        """Makes the partition of the dofs for prescribed displacements

        Makes the following partition of a given matrix::

            k = | kkk    kku |
                | kuk    kuu |

        Parameters
        ----------
        k : scipy.sparse.coo_matrix
            Matrix to be partitioned.
        return_kkk : bool, optional
            If the region `kkk` must be returned.
        return_kku : bool, optional
            If the region `kku` must be returned.
        return_kuk : bool, optional
            If the region `kuk` must be returned.

        Returns
        -------
        out : dict
            A ``dict`` object containing the keys for the
            corresponding sub-matrices ``kkk``, ``kku``, ``kuk``, ``kuu``.
            The sub-matrix ``out['kuu']`` is a ``scipy.sparse.csr_matrix``,
            while the others are 2-D ``np.ndarray`` objects.

        """
        if not isinstance(k, coo_matrix):
            k = coo_matrix(k)

        if return_kkk:
            kkk = coo_matrix(np.zeros((self.num0, self.num0)))
            ind = np.where(((k.row < self.num0) & (k.col < self.num0)))[0]
            kkk.row = np.take(k.row, ind)
            kkk.col = np.take(k.col, ind)
            kkk.data = np.take(k.data, ind)
            kkk = kkk.toarray()
            kkk = np.delete(kkk, self.excluded_dofs, axis=0)
            kkk = np.delete(kkk, self.excluded_dofs, axis=1)

        if return_kku:
            kku = coo_matrix(np.zeros((self.num0, k.shape[0])))
            ind = np.where(k.row < self.num0)[0]
            kku.row = np.take(k.row, ind)
            kku.col = np.take(k.col, ind)
            kku.data = np.take(k.data, ind)
            kku = kku.toarray()
            kku = np.delete(kku, self.excluded_dofs, axis=1)

        if return_kuk:
            kuk = coo_matrix(np.zeros((k.shape[0], self.num0)))
            ind = np.where(k.col < self.num0)[0]
            kuk.row = np.take(k.row, ind)
            kuk.col = np.take(k.col, ind)
            kuk.data = np.take(k.data, ind)
            kuk = kuk.toarray()
            kuk = np.delete(kuk, self.excluded_dofs, axis=0)

        rows = np.sort(self.excluded_dofs)[::-1]
        cols = np.sort(self.excluded_dofs)[::-1]

        kuu = k.copy()

        for r in rows:
            ind = np.where(kuu.row != r)[0]
            kuu.row[kuu.row > r] -= 1
            kuu.row = np.take(kuu.row, ind)
            kuu.col = np.take(kuu.col, ind)
            kuu.data = np.take(kuu.data, ind)
            kuu._shape = (kuu._shape[0]-1, kuu._shape[1])

        for c in cols:
            ind = np.where(kuu.col != c)[0]
            kuu.col[kuu.col > c] -= 1
            kuu.row = np.take(kuu.row, ind)
            kuu.col = np.take(kuu.col, ind)
            kuu.data = np.take(kuu.data, ind)
            kuu._shape = (kuu._shape[0], kuu._shape[1]-1)

        kuu = csr_matrix(kuu)

        out = {}
        out['kuu'] = kuu
        if return_kkk:
            out['kkk'] = kkk
        if return_kku:
            out['kku'] = kku
        if return_kuk:
            out['kuk'] = kuk

        return out


    def calc_full_c(self, cu, inc=1.):
        """Returns the full set of Ritz constants

        When prescribed displacements take place the matrices and the Ritz
        constants are partitioned like::

            k = | kkk    kku |
                | kuk    kuu |

        and the corresponding Ritz constants::

            c = | ck |
                | cu |

        This function adds the set of known Ritz constants (``ck``)
        to the set of unknown (``cu``) based on the prescribed displacements.

        Parameters
        ----------
        cu : np.ndarray
            The set of unknown Ritz constants
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants.

        Returns
        -------
        c : np.ndarray
            The full set of Ritz constants.

        """
        c = cu.copy()
        size = self.get_size()
        if c.shape[0] == size:
            for dof in self.excluded_dofs:
                c[dof] *= inc
            return c
        ordered = sorted(zip(self.excluded_dofs,
                             self.excluded_dofs_ck), key=lambda x:x[0])
        for dof, cai in ordered:
            c = np.insert(c, dof, inc*cai)
        return c


    def _default_field(self, xs, ts, gridx, gridt):
        if xs is None or ts is None:
            xs = linspace(0, self.L, gridx)
            ts = linspace(-pi, pi, gridt)
            xs, ts = np.meshgrid(xs, ts, copy=False)
        xs = np.atleast_1d(np.array(xs, dtype=DOUBLE))
        ts = np.atleast_1d(np.array(ts, dtype=DOUBLE))
        xshape = xs.shape
        tshape = ts.shape
        if xshape != tshape:
            raise ValueError('Arrays xs and ts must have the same shape')
        self.Xs = xs
        self.Ts = ts
        xs = xs.ravel()
        ts = ts.ravel()

        return xs, ts, xshape, tshape


    def _calc_linear_matrices(self, combined_load_case=None):
        self._rebuild()
        msg('Calculating linear matrices... ', level=2)

        fk0, fk0_cyl, fkG0, fkG0_cyl, k0edges = modelDB.get_linear_matrices(
                                                    self, combined_load_case)
        model = self.model
        alpharad = self.alpharad
        cosa = self.cosa
        r1 = self.r1
        r2 = self.r2
        L = self.L
        m1 = self.m1
        m2 = self.m2
        n2 = self.n2
        s = self.s
        laminaprops = self.laminaprops
        plyts = self.plyts
        stack = self.stack
        P = self.P
        T = self.T
        E11 = self.E11
        nu = self.nu
        h = self.h

        Fc = self.Nxxtop[0]*(2*pi*r2*cosa)

        lam = self.lam

        if stack != [] and self.F_reuse is None:
            lam = laminate.read_stack(stack, plyts=plyts,
                                             laminaprops=laminaprops)

        if 'clpt' in model:
            if self.F_reuse is not None:
                msg('')
                msg('Reusing F matrix...', level=2)
                F = self.F_reuse
            elif lam is not None:
                F = lam.ABD
            else:
                F = self.F

        elif 'fsdt' in model:
            if self.F_reuse is not None:
                msg('')
                msg('Reusing F matrix...', level=2)
                F = self.F_reuse
            elif lam is not None:
                F = lam.ABDE
                F[6:, 6:] *= self.K
            else:
                F = self.F

        if self.force_orthotropic_laminate:
            msg('')
            msg('Forcing orthotropic laminate...', level=2)
            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

        self.lam = lam
        self.F = F

        if self.is_cylinder:
            if 'iso_' in model:
                k0 = fk0_cyl(r2, L, E11, nu, h, m1, m2, n2)
            else:
                k0 = fk0_cyl(r2, L, F, m1, m2, n2)

            if not combined_load_case:
                kG0 = fkG0_cyl(Fc, P, T, r2, L, m1, m2, n2)
            else:
                kG0_Fc = fkG0_cyl(Fc, 0, 0, r2, L, m1, m2, n2)
                kG0_P = fkG0_cyl(0, P, 0, r2, L, m1, m2, n2)
                kG0_T = fkG0_cyl(0, 0, T, r2, L, m1, m2, n2)
        else:
            if 'iso_' in model:
                k0 = fk0(alpharad, r2, L, E11, nu, h, m1, m2, n2, s)
            else:
                k0 = fk0(alpharad, r2, L, F, m1, m2, n2, s)
            if not combined_load_case:
                kG0 = fkG0(Fc, P, T, r2, alpharad, L, m1, m2, n2, s)
            else:
                kG0_Fc = fkG0(Fc, 0, 0, r2, alpharad, L, m1, m2, n2, s)
                kG0_P = fkG0(0, P, 0, r2, alpharad, L, m1, m2, n2, s)
                kG0_T = fkG0(0, 0, T, r2, alpharad, L, m1, m2, n2, s)

        if k0edges is not None:
            k0 = csr_matrix(k0) + csr_matrix(k0edges)

        assert np.any((np.isnan(k0.data) | np.isinf(k0.data))) == False

        k0 = make_symmetric(k0)

        if not combined_load_case:
            assert np.any((np.isnan(kG0.data) | np.isinf(kG0.data))) == False
            self.kG0 = make_symmetric(kG0)
        else:
            assert np.any((np.isnan(kG0_Fc.data)
                           | np.isinf(kG0_Fc.data))) == False
            assert np.any((np.isnan(kG0_P.data)
                           | np.isinf(kG0_P.data))) == False
            assert np.any((np.isnan(kG0_T.data)
                           | np.isinf(kG0_T.data))) == False
            self.kG0_Fc = make_symmetric(kG0_Fc)
            self.kG0_P = make_symmetric(kG0_P)
            self.kG0_T = make_symmetric(kG0_T)

        k = self.exclude_dofs_matrix(k0, return_kuk=True)
        k0uk = k['kuk']
        k0uu = k['kuu']

        self.k0 = k0
        self.k0uk = k0uk
        self.k0uu = k0uu

        #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)


    def calc_k0(self):
        if self.k0uu is None:
            self._calc_linear_matrices()
        return self.k0uu


    def calc_kT(self, c, inc=1.):
        r"""Calculates the tangent stiffness matrix

        The following attributes will affect the numerical integration:

        =================    ================================================
        Attribute            Description
        =================    ================================================
        ``ni_num_cores``     ``int``, number of cores used for the numerical
                             integration
        ``ni_method``        ``str``, integration method:
                                 - ``'trapz2d'`` for 2-D Trapezoidal's rule
                                 - ``'simps2d'`` for 2-D Simpsons' rule
        ``nx``               ``int``, number of integration points along the
                             `x` coordinate
        ``nt``               ``int``, number of integration points along the
                             `\theta` coordinate
        =================    ================================================

        Parameters
        ----------
        c : np.ndarray
            The Ritz constants vector of the current state.
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants using :meth:`calc_full_c`.

        Returns
        -------
        kTuu : sparse matrix
            The tangent stiffness matrix corresponding to the unknown degrees
            of freedom.

        """
        self._calc_NL_matrices(c, inc=inc)
        return self.kTuu


    def lb(self, c=None, tol=0, combined_load_case=None):
        """Performs a linear buckling analysis

        The following parameters of the ``ConeCyl`` object will affect the
        linear buckling analysis:

        =======================    =====================================
        Attribute                  Description
        =======================    =====================================
        ``num_eigenvalues``        Number of eigenvalues to be extracted
        ``num_eigvalues_print``    Number of eigenvalues to print after
                                   the analysis is completed
        =======================    =====================================

        Parameters
        ----------
        combined_load_case : int, optional
            It tells whether the linear buckling analysis must be computed
            considering combined load cases, each value will tell
            the algorithm to rearrange the linear matrices in a different
            way. The valid values are ``1``, or ``2``, where:

            - ``1`` : find the critical axial load for a fixed torsion load
            - ``2`` : find the critical axial load for a fixed pressure load
            - ``3`` : find the critical torsion load for a fixed axial load

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

        """
        model_dict = get_model(self.model)
        if not model_dict['linear buckling']:
            msg('________________________________________________')
            msg('')
            warn('Model {} cannot be used in linear buckling analysis!'.
                 format(self.model))
            msg('________________________________________________')

        msg('Running linear buckling analysis...')

        if self.Fc is None and self.Nxxtop is None:
            warn('using Fc = 1.', level=1)
            self.Fc = 1.
        if self.pdC is None:
            self.pdC = False

        self._calc_linear_matrices(combined_load_case=combined_load_case)

        #TODO maybe a better estimator to sigma would be to run
        #     a preliminary eigsh using a small m2 and n2
        #NOTE runs faster for self.k0 than -self.k0, so that the negative
        #     sign is applied later

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

        model_dict = get_model(self.model)
        i0 = model_dict['i0']
        num0 = model_dict['num0']
        num1 = model_dict['num1']
        num2 = model_dict['num2']

        pos = num0

        if not combined_load_case:
            M = csr_matrix(self.k0)
            A = csr_matrix(self.kG0)
        elif combined_load_case == 1:
            M = csr_matrix(self.k0) + csr_matrix(self.kG0_T)
            A = csr_matrix(self.kG0_Fc)
        elif combined_load_case == 2:
            M = csr_matrix(self.k0) + csr_matrix(self.kG0_P)
            A = csr_matrix(self.kG0_Fc)
        elif combined_load_case == 3:
            M = csr_matrix(self.k0) + csr_matrix(self.kG0_Fc)
            A = csr_matrix(self.kG0_T)

        A = A[pos:, pos:]
        M = M[pos:, pos:]

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

        eigvals = -1./eigvals
        self.eigvals = eigvals
        self.eigvecs = np.vstack((np.zeros((pos, self.num_eigvalues)),
                                  eigvecs))

        msg('finished!', level=2)

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


    def eigen(self, c=None, tol=0, kL=None, kG=None):
        """Performs a non-linear eigenvalue analysis at a given state

        The following attributes of the ``ConeCyl`` object will affect the
        non-linear eigenvalue analysis:

        =======================    =====================================
        Attribute                  Description
        =======================    =====================================
        ``num_eigenvalues``        Number of eigenvalues to be extracted
        ``num_eigvalues_print``    Number of eigenvalues to print after
                                   the analysis is completed
        =======================    =====================================

        Additionally, the non-linear analysis parameters described in
        :meth:`static` will affect the integration of the non-linear matrices
        ``kL`` and ``kG`` if they are not given as input parameters.

        Parameters
        ----------
        combined_load_case : int, optional
            It tells whether the linear buckling analysis must be computed
            considering combined load cases, each value will tell
            the algorithm to rearrange the linear matrices in a different
            way. The valid values are ``1``, or ``2``, where:

            - ``1`` : find the critical axial load for a fixed torsion load
            - ``2`` : find the critical axial load for a fixed pressure load
            - ``3`` : find the critical torsion load for a fixed axial load

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

        """
        model_dict = get_model(self.model)
        if not model_dict['linear buckling']:
            msg('________________________________________________')
            msg('')
            warn('Model {} cannot be used in linear buckling analysis!'.
                 format(self.model))
            msg('________________________________________________')

        msg('Running linear buckling analysis...')

        if self.Fc is None:
            self.Fc = 1.
        if self.pdC is None:
            self.pdC = False

        self._calc_linear_matrices(combined_load_case=combined_load_case)

        #TODO maybe a better estimator to sigma would be to run
        #     a preliminary eigsh using a small m2 and n2
        #NOTE runs faster for self.k0 than -self.k0, so that the negative
        #     sign is applied later
        msg('Eigenvalue solver... ', level=2)

        model_dict = get_model(self.model)
        i0 = model_dict['i0']
        num0 = model_dict['num0']
        num1 = model_dict['num1']
        num2 = model_dict['num2']

        pos = num0

        if not combined_load_case:
            M = csr_matrix(self.k0)
            A = csr_matrix(self.kG0)
        elif combined_load_case == 1:
            M = csr_matrix(self.k0) + csr_matrix(self.kG0_T)
            A = csr_matrix(self.kG0_Fc)
        elif combined_load_case == 2:
            M = csr_matrix(self.k0) + csr_matrix(self.kG0_P)
            A = csr_matrix(self.kG0_Fc)
        elif combined_load_case == 3:
            M = csr_matrix(self.k0) + csr_matrix(self.kG0_Fc)
            A = csr_matrix(self.kG0_T)

        A = A[pos:, pos:]
        M = M[pos:, pos:]

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

        eigvals = (-1./eigvals)
        self.eigvals = eigvals
        self.eigvecs = np.vstack((np.zeros((pos, self.num_eigvalues)),
                                  eigvecs))

        msg('finished!', level=2)

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


    def _calc_NL_matrices(self, c, inc=1., with_kLL=None, with_k0L=None):
        r"""Calculates the non-linear stiffness matrices

        Parameters
        ----------
        c : np.ndarray
            Ritz constants representing the current state to calculate the
            stiffness matrices.
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants using :meth:`calc_full_c`.
        with_kLL : bool, optional
            When ``with_kLL=False`` assumes kLL << than k0L and kG.
        with_k0L : bool, optional
            When ``with_k0L=False`` assumes k0L << than kLL and kG.

        Notes
        -----
        Nothing is returned, the calculated matrices

        """
        c = self.calc_full_c(c, inc=inc)

        if self.k0 is None:
            self._calc_linear_matrices()
        if with_k0L is None:
            with_k0L = self.with_k0L
        if with_kLL is None:
            with_kLL = self.with_kLL

        msg('Calculating non-linear matrices...', level=2)
        alpharad = self.alpharad
        r2 = self.r2
        L = self.L
        tLArad = self.tLArad
        F = self.F
        m1 = self.m1
        m2 = self.m2
        n2 = self.n2
        c0 = self.c0
        m0 = self.m0
        n0 = self.n0
        funcnum = self.funcnum

        model = self.model
        model_dict = get_model(model)

        nlmodule = model_dict['non-linear']
        ni_method = self.ni_method
        num_cores = self.ni_num_cores
        nx = self.nx
        nt = self.nt

        if nlmodule:
            calc_k0L = nlmodule.calc_k0L
            calc_kLL = nlmodule.calc_kLL
            if 'iso_' in model:
                calc_kG = modelDB.db[model[4:]]['non-linear'].calc_kG
            else:
                calc_kG = nlmodule.calc_kG

            kG = calc_kG(c, alpharad, r2, L, tLArad, F, m1, m2, n2, nx=nx,
                    nt=nt, num_cores=num_cores, method=ni_method, c0=c0,
                    m0=m0, n0=n0)
            kG = make_symmetric(kG)

            if 'iso_' in model:
                E11 = self.E11
                nu = self.nu
                h = self.h
                if with_k0L:
                    k0L = calc_k0L(c, alpharad, r2, L, tLArad, E11, nu, h, m1,
                            m2, n2, nx=nx, nt=nt, num_cores=num_cores,
                            method=ni_method, c0=c0, m0=m0, n0=n0)
                else:
                    k0L = kG*0
                if with_kLL:
                    kLL = calc_kLL(c, alpharad, r2, L, tLArad, E11, nu, h, m1,
                            m2, n2, nx=nx, nt=nt, num_cores=num_cores,
                            method=ni_method, c0=c0, m0=m0, n0=n0)
                    kLL = make_symmetric(kLL)

                else:
                    kLL = kG*0

            else:
                if with_k0L:
                    k0L = calc_k0L(c, alpharad, r2, L, tLArad, F, m1, m2, n2,
                            nx=nx, nt=nt, num_cores=num_cores,
                            method=ni_method, c0=c0, m0=m0, n0=n0)
                else:
                    k0L = kG*0
                if with_kLL:
                    kLL = calc_kLL(c, alpharad, r2, L, tLArad, F, m1, m2, n2,
                            nx=nx, nt=nt, num_cores=num_cores,
                            method=ni_method, c0=c0, m0=m0, n0=n0)
                    kLL = make_symmetric(kLL)

                else:
                    kLL = kG*0

        else:
            raise ValueError(
            'Non-Linear analysis not implemented for model {0}'.format(model))

        kL0 = k0L.T

        #TODO maybe slow...
        kT = coo_matrix(self.k0 + k0L + kL0 + kLL + kG)
        # kS was deprecated, now fint is integrated numerically
        #kS = coo_matrix(self.k0 + k0L/2 + kL0 + kLL/2)

        k = self.exclude_dofs_matrix(kT, return_kuk=True)
        self.kTuk = k['kuk']
        self.kTuu = k['kuu']

        #NOTE intended for non-linear eigenvalue analyses
        self.kL = csr_matrix(self.k0 + k0L + kL0 + kLL)
        self.kG = csr_matrix(kG)

        msg('finished!', level=2)


    def uvw(self, c, xs=None, ts=None, gridx=300, gridt=300, inc=1.):
        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``, ``phit`` of the ``ConeCyl`` 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.
        ts : np.ndarray
            The ``theta`` 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.
        gridt : int
            Number of points along the `theta` where to calculate the
            displacement field.
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants using :meth:`calc_full_c`.

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

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

        """
        xs, ts, xshape, tshape = self._default_field(xs, ts, gridx, gridt)
        alpharad = self.alpharad
        tLArad = self.tLArad
        m1 = self.m1
        m2 = self.m2
        n2 = self.n2
        r2 = self.r2
        L = self.L

        c = self.calc_full_c(c, inc=inc)

        model_dict = get_model(self.model)
        fuvw = model_dict['commons'].fuvw
        us, vs, ws, phixs, phits = fuvw(c, m1, m2, n2, alpharad, r2, L,
                                        tLArad, xs, ts, 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.phit = phits.reshape(xshape)

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


    def strain(self, c, xs=None, ts=None, gridx=300, gridt=300, inc=1.):
        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.
        ts : np.ndarray, optional
            The `\theta` coordinates where to calculate the strains, must
            have the same shape as ``xs``.
        gridx : int, optional
            When ``xs`` and ``ts`` are not supplied, ``gridx`` and ``gridt``
            are used.
        gridt : int, optional
            When ``xs`` and ``ts`` are not supplied, ``gridx`` and ``gridt``
            are used.
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants using :meth:`calc_full_c`.

        """
        xs, ts, xshape, tshape = self._default_field(xs, ts, gridx, gridt)

        alpharad = self.alpharad
        L = self.L
        r2 = self.r2
        sina = self.sina
        cosa = self.cosa
        tLArad = self.tLArad
        m1 = self.m1
        m2 = self.m2
        n2 = self.n2
        c0 = self.c0
        m0 = self.m0
        n0 = self.n0
        funcnum = self.funcnum
        model = self.model
        model_dict = get_model(model)
        NL_kinematics = model.split('_')[1]
        fstrain = model_dict['commons'].fstrain
        e_num = model_dict['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))

        c = self.calc_full_c(c, inc=inc)

        es = fstrain(c, sina, cosa, tLArad, xs, ts, r2, L,
                     m1, m2, n2, c0, m0, n0, funcnum, int_NL_kinematics,
                     self.out_num_cores)

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


    def stress(self, c, xs=None, ts=None, gridx=300, gridt=300, inc=1.):
        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.
        ts : np.ndarray, optional
            The `\theta` coordinates where to calculate the strains, must
            have the same shape as ``xs``.
        gridx : int, optional
            When ``xs`` and ``ts`` are not supplied, ``gridx`` and ``gridt``
            are used.
        gridt : int, optional
            When ``xs`` and ``ts`` are not supplied, ``gridx`` and ``gridt``
            are used.
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants using :meth:`calc_full_c`.

        """
        xs, ts, xshape, tshape = self._default_field(xs, ts, gridx, gridt)

        F = self.F
        alpharad = self.alpharad
        L = self.L
        r2 = self.r2
        sina = self.sina
        cosa = self.cosa
        tLArad = self.tLArad
        m1 = self.m1
        m2 = self.m2
        n2 = self.n2
        c0 = self.c0
        m0 = self.m0
        n0 = self.n0
        funcnum = self.funcnum
        model = self.model
        model_dict = get_model(model)
        NL_kinematics = model.split('_')[1]
        fstress = model_dict['commons'].fstress
        e_num = model_dict['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))

        c = self.calc_full_c(c, inc=inc)

        Ns = fstress(c, F, sina, cosa, tLArad, xs, ts, r2, L,
                     m1, m2, n2, c0, m0, n0, funcnum, int_NL_kinematics,
                     self.out_num_cores)
        return Ns.reshape((xshape + (e_num,)))


    def calc_fint(self, c, inc=1., m=1, return_u=True):
        r"""Calculates the internal force vector `\{F_{int}\}`

        The following attributes will affect the numerical integration:

        =================    ================================================
        Attribute            Description
        =================    ================================================
        ``ni_num_cores``     ``int``, number of cores used for the numerical
                             integration
        ``ni_method``        ``str``, integration method:
                                 - ``'trapz2d'`` for 2-D Trapezoidal's rule
                                 - ``'simps2d'`` for 2-D Simpsons' rule
        ``nx``               ``int``, number of integration points along the
                             `x` coordinate
        ``nt``               ``int``, number of integration points along the
                             `\theta` coordinate
        =================    ================================================

        Parameters
        ----------
        c : np.ndarray
            The Ritz constants that will be used to compute the internal
            forces.
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants using :meth:`calc_full_c`.
        m : integer, optional
            A multiplier to be applied to ``nx`` and ``nt``, if one
            whishes to use more integration points.
        return_u : bool, optional
            If the internal force vector corresponsing to the unknown
            set of Ritz constants should be returned.

        Returns
        -------
        fint : np.ndarray
            The internal force vector.

        """
        c = self.calc_full_c(c, inc=inc)
        if 'iso_' in self.model:
            nlmodule = modelDB.db[self.model[4:]]['non-linear']
        else:
            nlmodule = modelDB.db[self.model]['non-linear']
        ni_method = self.ni_method
        ni_num_cores = self.ni_num_cores
        nx = self.nx*m
        nt = self.nt*m
        fint = nlmodule.calc_fint_0L_L0_LL(c, self.alpharad, self.r2, self.L,
                self.tLArad, self.F, self.m1, self.m2, self.n2, nx, nt,
                ni_num_cores, ni_method, self.c0, self.m0, self.n0)
        fint += self.k0*c
        if return_u:
            fint = np.delete(fint, self.excluded_dofs)

        return fint


    def add_SPL(self, PL, pt=0.5, thetadeg=0., increment=False):
        """Add a Single Perturbation Load `\{{F_{PL}}_i\}`

        Adds a perturbation load to the ``ConeCyl`` object, the perturbation
        load is a particular case of the punctual load with only a normal
        component.

        Parameters
        ----------
        PL : float
            The perturbation load value.
        pt : float, optional
            The normalized position along the `x` axis in which the new SPL
            will be included.
        thetadeg : float, optional
            The angular position of the SPL in degrees.
        increment : bool, optional
            If this perturbation load should be incrementally applied in a
            non-linear analysis.

        Notes
        -----
        Each single perturbation load is added to the ``forces`` parameter
        of the ``ConeCyl`` object, which may be changed by the analyst at
        any time.

        """
        self._rebuild()
        thetarad = deg2rad(thetadeg)
        if increment:
            self.forces_inc.append([pt*self.L, thetarad, 0., 0., PL])
        else:
            self.forces.append([pt*self.L, thetarad, 0., 0., PL])


    def add_force(self, x, thetadeg, fx, ftheta, fz, increment=False):
        r"""Add a punctual force

        Adds a force vector `\{f_x, f_\theta, f_z\}^T` to the ``forces``
        parameter of the ``ConeCyl`` object.

        Parameters
        ----------
        x : float
            The `x` position.
        thetadeg : float
            The `\theta` position in degrees.
        fx : float
            The `x` component of the force vector.
        ftheta : float
            The `\theta` component of the force vector.
        fz : float
            The `z` component of the force vector.
        increment : bool, optional
            If this punctual force should be incrementally applied in a
            non-linear analysis.

        """
        thetarad = deg2rad(thetadeg)
        if incremented:
            self.forces_inc.append([x, thetarad, fx, ftheta, fz])
        else:
            self.forces.append([x, thetarad, fx, ftheta, fz])


    def calc_fext(self, inc=1., kuk=None, 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\}`.

        kuk : np.ndarray, optional
            Obsolete, created for displacement controlled analyses, but the
            implementation has not been finished, see
            :meth:`exclude_dofs_matrix`.
        silent : bool, optional
            A boolean to tell whether the msg messages should be printed.

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

        """
        self._rebuild()
        if self.k0 is None:
            self._calc_linear_matrices()

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


        uTM = inc*self.uTM
        Nxxtop = inc*self.Nxxtop
        thetaTrad = inc*self.thetaTrad

        sina = self.sina
        cosa = self.cosa
        r2 = self.r2
        L = self.L
        tLArad = self.tLArad
        m1 = self.m1
        m2 = self.m2
        n2 = self.n2
        pdC = self.pdC
        pdT = self.pdT
        model = self.model

        model_dict = get_model(model)

        i0 = model_dict['i0']
        j0 = model_dict['j0']
        num0 = model_dict['num0']
        num1 = model_dict['num1']
        num2 = model_dict['num2']
        dofs = model_dict['dofs']
        fg = model_dict['commons'].fg

        size = self.get_size()

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

        fext = np.delete(fext, self.excluded_dofs)

        # constant punctual forces
        for i, force in enumerate(self.forces):
            x, theta, fx, ftheta, fz = force
            fg(g, m1, m2, n2, r2, x, theta, L, cosa, tLArad)

            gu = np.delete(g, self.excluded_dofs, axis=1)

            if dofs == 3:
                fpt = np.array([[fx, ftheta, fz]])
            elif dofs == 5:
                fpt = np.array([[fx, ftheta, fz, 0, 0]])
            fext += -fpt.dot(gu).ravel()

        # incremented punctual forces
        for i, force in enumerate(self.forces_inc):
            x, theta, fx, ftheta, fz = force
            fg(g, m1, m2, n2, r2, x, theta, L, cosa, tLArad)

            gu = np.delete(g, self.excluded_dofs, axis=1)

            if dofs == 3:
                fpt = inc*np.array([[fx, ftheta, fz]])
            elif dofs == 5:
                fpt = inc*np.array([[fx, ftheta, fz, 0, 0]])
            fext += -fpt.dot(gu).ravel()

        # axial load
        fext_tmp = np.zeros(size, dtype=DOUBLE)
        if not 0 in self.excluded_dofs:
            fext_tmp[0] += Nxxtop[0]*(2*pi*r2)/cosa
            if 'bc2' in model or 'bc4' in model:
                for j2 in range(j0, n2+j0):
                    for i2 in range(i0, m2+i0):
                        row = (num0 + num1*m1
                               + (i2-i0)*num2 + (j2-j0)*num2*m2)
                        rowNxx = 1+2*(j2-j0)
                        fext_tmp[row+0]+=(Nxxtop[rowNxx+0]*pi*r2)
                        fext_tmp[row+1]+=(Nxxtop[rowNxx+1]*pi*r2)
        else:
            if kuk is None:
                kuk_C = self.k0uk[:, 0].ravel()
            else:
                kuk_C = kuk[:, 0].ravel()
            fext += -uTM*kuk_C

        if not 2 in self.excluded_dofs:
            fext_tmp[2] += Nxxtop[2]*(2*pi*r2)/cosa

        # pressure
        P = self.P + inc*self.P_inc
        if P != 0:
            if 'clpt' in model:
                for i1 in range(i0, m1+i0):
                    if i1 == 0:
                        continue
                    col = num0 + (i1-i0)*num1
                    fext_tmp[col+2] += P*(L*2./i1*(r2 - (-1)**i1*(r2 + L*sina)))
            elif 'fsdt' in model:
                #TODO it might be the same as for the CLPT
                raise NotImplementedError(
                    'Pressure not implemented for static analysis for FSDT')

        fext_tmp = np.delete(fext_tmp, self.excluded_dofs)
        fext += fext_tmp

        # torsion
        if pdT:
            if kuk is None:
                kuk_T = self.k0uk[:, 1].ravel()
            else:
                kuk_T = kuk[:, 1].ravel()
            fext += -thetaTrad*kuk_T
        else:
            T = self.T + inc*self.T_inc
            if T != 0:
                fg(g, m1, m2, n2, r2, 0, 0, L, cosa, tLArad)
                gu = np.delete(g, self.excluded_dofs, axis=1)
                if dofs == 3:
                    fpt = np.array([[0, T/r2, 0]])
                elif dofs == 5:
                    fpt = np.array([[0, T/r2, 0, 0, 0]])
                fext += fpt.dot(gu).ravel()

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

        return fext


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

        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 msg 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 the ``cs`` parameter of the
        ``ConeCyl`` object. The actual increments used in the non-linear
        analysis are stored in the ``increments`` parameter.

        """
        self.cs = []
        self.increments = []

        if self.pdC:
            text ='Non-linear analysis with prescribed displacements'
            raise NotImplementedError(text)

        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]['linear static']:
            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)
            lob('________________________________________________',
                level=1, silent=silent)
            raise

        self.analysis.static(NLgeom=NLgeom, silent=silent)
        self.cs = self.analysis.cs
        self.increments = self.analysis.increments

        return self.cs


    def plot(self, c, invert_x=False, plot_type=1, vec='w',
             deform_u=False, deform_u_sf=100.,
             filename='',
             ax=None, figsize=(3.5, 2.), save=True,
             add_title=True, title='',
             colorbar=False, cbar_nticks=2, cbar_format=None,
             cbar_title='', cbar_fontsize=10,
             aspect='equal', clean=True, dpi=400,
             texts=[], xs=None, ts=None, gridx=300, gridt=300,
             num_levels=400, inc=1.):
        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'``, ``'phit'``,
                            ``'magnitude'``
            - Strain: ``'exx'``, ``'ett'``, ``'gxt'``, ``'kxx'``, ``'ktt'``,
              ``'kxt'``, ``'gtz'``, ``'gxz'``
            - Stress: ``'Nxx'``, ``'Ntt'``, ``'Nxt'``, ``'Mxx'``, ``'Mtt'``,
              ``'Mxt'``, ``'Qt'``, ``'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_x : bool, optional
            Inverts the `x` axis of the plot. It may be used to match
            the coordinate system of the finite element models created
            using the ``desicos.abaqus`` module.
        plot_type : int, optional
            For cylinders only ``4`` and ``5`` are valid.
            For cones all the following types can be used:

            - ``1``: concave up (with ``invert_x=False``) (default)
            - ``2``: concave down (with ``invert_x=False``)
            - ``3``: stretched closed
            - ``4``: stretched opened (`r \times \theta` vs. `H`)
            - ``5``: stretched opened (`\theta` vs. `H`)

        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 ``ConeCyl`` 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)``.
        add_title : bool, optional
            If a title should be added to the figure.
        title : str, optional
            If any string is given ``add_title`` will be ignored and the given
            title added 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.
        ts : np.ndarray, optional
            The ``theta`` 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.
        gridt : int, optional
            Number of points along the `theta` where to calculate the
            displacement field.
        num_levels : int, optional
            Number of contour levels (higher values make the contour smoother).
        inc : float, optional
            Load increment, necessary to calculate the full set of Ritz
            constants using :meth:`calc_full_c`.

        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, phitbkp = (self.u, self.v, self.w,
                                              self.phix, self.phit)

        import matplotlib.pyplot as plt
        import matplotlib

        from plotutils import get_filename

        c = self.calc_full_c(c, inc=inc)

        msg('Computing field variables...', level=1)
        displs = ['u', 'v', 'w', 'phix', 'phit', 'magnitude', 'test']
        strains = ['exx', 'ett', 'gxt', 'kxx', 'ktt', 'kxt', 'gtz', 'gxz']
        stresses = ['Nxx', 'Ntt', 'Nxt', 'Mxx', 'Mtt', 'Mxt', 'Qt', 'Qx']
        if vec in displs or 'eq_' in vec:
            self.uvw(c, xs=xs, ts=ts, gridx=gridx, gridt=gridt, inc=inc)
            if vec == 'magnitude':
                u = self.u
                v = self.v
                w = self.w
                field = (u**2 + v**2 + w**2)**0.5
            elif 'eq_' in vec:
                u = self.u
                v = self.v
                w = self.w
                field = eval(vec[3:])
            else:
                field = getattr(self, vec)
        elif vec in strains:
            es = self.strain(c, xs=xs, ts=ts,
                             gridx=gridx, gridt=gridt, inc=inc)
            field = es[..., strains.index(vec)]
        elif vec in stresses:
            Ns = self.stress(c, xs=xs, ts=ts,
                             gridx=gridx, gridt=gridt, inc=inc)
            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
        Ts = self.Ts

        vecmin = field.min()
        vecmax = field.max()

        levels = np.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')

        def r(x):
            return self.r2 + x*self.sina

        if self.is_cylinder:
            plot_type=4
        if plot_type == 1:
            r_plot = self.r2/self.sina + Xs
            r_plot_max = self.r2/self.sina + self.L
            y = r_plot_max - r_plot*cos(Ts*self.sina)
            x = r_plot*sin(Ts*self.sina)
        elif plot_type == 2:
            r_plot = self.r2/self.sina + Xs
            y = r_plot*cos(Ts*self.sina)
            x = r_plot*sin(Ts*self.sina)
        elif plot_type == 3:
            r_plot = self.r2/self.sina + Xs
            r_plot_max = self.r2/self.sina + self.L
            y = r_plot_max - r_plot*cos(Ts)
            x = r_plot*sin(Ts)
        elif plot_type == 4:
            x = r(Xs)*Ts
            y = self.L-Xs
        elif plot_type == 5:
            x = Ts
            y = Xs
        if deform_u:
            if vec in displs:
                pass
            else:
                self.uvw(c, xs=xs, ts=ts, gridx=gridx, gridt=gridt, inc=inc)
            field_u = self.u
            y -= deform_u_sf*field_u
        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 = np.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_x:
            ax.invert_yaxis()

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

        elif add_title:
            if self.analysis.last_analysis == 'static':
                ax.set_title(r'$m_1={0}$, $m_2={1}$, $n_2={2}$'.
                             format(self.m1, self.m2, self.n2))

            elif self.analysis.last_analysis == 'lb':
                ax.set_title(
           r'$m_1={0}$, $m2={1}$, $n2={2}$, $\lambda_{{CR}}={3:1.3e}$'.format(
                self.m1, self.m2, self.n2, self.eigvals[0]))

        fig.tight_layout()
        ax.set_aspect(aspect)

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

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

        if save:
            if not filename:
                filename = get_filename(self)
            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 phitbkp is not None:
            self.phit = phitbkp

        msg('finished!')

        return ax


    def plotAbaqus(self, frame, fieldOutputKey, vec, nodes, numel_cir,
                   elem_type='S4R', ignore=[],
                   ax=None, figsize=(3.3, 3.3), save=True,
                   aspect='equal', clean=True, plot_type=1,
                   outpath='', filename='', npzname='', pyname='',
                   num_levels=400):
        r"""Print a field output for a cylinder/cone model from Abaqus

        This function is intended to be used with models created by the
        `DESICOS plug-in for Abaqus <http://desicos.github.io/desicos/>`_,
        where a **mapped mesh** is used and the models comparable to the
        models of :mod:`compmech.conecyl`.

        The ``frame`` and ``nodes`` input types are described in
        `Abaqus Scripting Reference Manual
        <http://abaqus.me.chalmers.se/v6.11/books/ker/default.htm>`_ and
        they can be obtained through:

        >>> frame = session.odbs['odb_name.odb'].steps['step_name'].frames[0]
        >>> nodes = mdb.models['model_name'].parts['part_name'].nodes

        Parameters
        ----------
        frame : OdbFrame
            The frame from where the field output will be taken from.
        fieldOutputKey : str
            The field output key to be used. It must be available in
            ``frame.fieldOutputs.keys()``. This function was tested with
            ``'UT'`` and ``'U'`` only.
        vec : str
            The displacement vector to be plotted:
            ``'u'``, ``'v'`` or ``'w'``.
        nodes : MeshNodeArray
            The part nodes.
        numel_cir : int
            The number of elements around the circumference.
        elem_type : str, optional
            The element type. The elements ``'S4R', 'S4R5'`` where tested.
        ignore : list, optional
            A list with the node ids to be ignored. It must contain any nodes
            outside the mapped mesh included in ``parts['part_name'].nodes``.
        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)``.
        save : bool, optional
            Flag telling whether the contour should be saved to an image file.
        aspect : str, optional
            String that will be passed to the ``AxesSubplot.set_aspect()``
            method.
        clean : bool, optional
            Clean axes ticks, grids, spines etc.
        plot_type : int, optional
            See :meth:`plot`.
        outpath : str, optional
            Output path where the data from Abaqus and the plots are
            saved (see notes).
        filename : str, optional
            The file name for the generated image file.
        npzname : str, optional
            The file name for the generated npz file.
        pyname : str, optional
            The file name for the generated Python file.
        num_levels : int, optional
            Number of contour levels (higher values make the contour smoother).

        Returns
        -------
        out : tuple
            Where ``out[0]`` and ``out[1]`` contain the circumferential and
            meridional grids of coordinates and ``out[2]`` the corresponding
            field output.

        Notes
        -----
        The data is saved using ``np.savez()`` into ``outpath`` as
        ``abaqus_output.npz`` with an accompanying script for plotting
        ``abaqus_output_plot.py``, very handy when Matplotlib is not
        importable from Abaqus.

        """
        workingplt = True
        if not npzname:
            npzname = 'abaqus_output.npz'
        npzname = os.path.join(outpath, npzname)
        if not pyname:
            pyname = 'abaqus_output_plot.py'
        pyname = os.path.join(outpath, pyname)
        if not filename:
            filename = 'plot_from_abaqus.png'
        filename = os.path.join(outpath, filename)
        try:
            import matplotlib.pyplot as plt
            import matplotlib
        except:
            workingplt = False
        try:
            if not frame:
                frame = utils.get_current_frame()
            if not frame:
                raise ValueError('A frame must be selected!')
            frame_num = int(frame.frameValue)
            coords = np.array([n.coordinates for n in nodes
                               if n.label not in ignore])
            #TODO include more outputs like stress etc
            field = frame.fieldOutputs[fieldOutputKey]

            uvw_rec = np.array([val.data for val in field.values
                if getattr(val.instance, 'name', None) == 'INSTANCECYLINDER'])
            u_rec = uvw_rec[:,0]
            v_rec = uvw_rec[:,1]
            w_rec = uvw_rec[:,2]

            res_alpha = np.arctan2(v_rec, u_rec)

            thetas = np.arctan2(coords[:, 1], coords[:, 0])

            sina = sin(self.alpharad)
            cosa = cos(self.alpharad)

            ucyl = -w_rec
            vcyl = v_rec*cos(thetas) - u_rec*sin(thetas)
            wcyl = v_rec*sin(thetas) + u_rec*cos(thetas)
            u = wcyl*sina + ucyl*cosa
            v = vcyl
            w = wcyl*cosa - ucyl*sina

            displ_vecs = {'u':u, 'v':v, 'w':w}
            uvw = displ_vecs[vec]

            zs = coords[:, 2]

            nt = numel_cir
            if 'S8' in elem_type:
                raise NotImplementedError('{0} elements!'.format(elem_type))
                #nt *= 2

            #first sort
            asort = zs.argsort()
            zs = zs[asort].reshape(-1, nt)
            thetas = thetas[asort].reshape(-1, nt)
            uvw = uvw[asort].reshape(-1, nt)

            #second sort
            asort = thetas.argsort(axis=1)
            for i, asorti in enumerate(asort):
                zs[i,:] = zs[i,:][asorti]
                thetas[i,:] = thetas[i,:][asorti]
                uvw[i,:] = uvw[i,:][asorti]

            H = self.H
            r2 = self.r2
            r1 = self.r1
            L = H/cosa

            def fr(z):
                return r1 - z*sina/cosa

            if self.alpharad == 0.:
                plot_type=4
            if plot_type == 1:
                r_plot = fr(zs)
                if self.alpharad == 0.:
                    r_plot_max = L
                else:
                    r_plot_max = r2/sina + L
                y = r_plot_max - r_plot*cos(thetas*sina)
                x = r_plot*sin(thetas*sina)
            elif plot_type == 2:
                r_plot = fr(zs)
                y = r_plot*cos(thetas*sina)
                x = r_plot*sin(thetas*sina)
            elif plot_type == 3:
                r_plot = fr(zs)
                r_plot_max = r2/sina + L
                y = r_plot_max - r_plot*cos(thetas)
                x = r_plot*sin(thetas)
            elif plot_type == 4:
                x = fr(zs)*thetas
                y = zs
            elif plot_type == 5:
                x = thetas
                y = zs

            cir = x
            mer = y
            field = uvw

            if workingplt:
                levels = np.linspace(field.min(), field.max(), 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')
                ax.contourf(cir, mer, field, levels=levels)
                ax.grid(False)
                ax.set_aspect(aspect)
                ax.xaxis.set_ticks_position('bottom')
                ax.yaxis.set_ticks_position('left')
                #lim = self.r2*pi
                #ax.xaxis.set_ticks([-lim, 0, lim])
                #ax.xaxis.set_ticklabels([r'$-\pi$', '$0$', r'$+\pi$'])
                #ax.set_title(
                    #r'$PL=20 N$, $F_{{C}}=50 kN$, $w_{{PL}}=\beta$, $mm$')
                if clean:
                    ax.xaxis.set_ticks_position('none')
                    ax.yaxis.set_ticks_position('none')
                    ax.xaxis.set_ticklabels([])
                    ax.yaxis.set_ticklabels([])
                    ax.set_frame_on(False)
                if save:
                    msg('Plot saved at: {0}'.format(filename))
                    plt.tight_layout()
                    plt.savefig(filename, transparent=True,
                                bbox_inches='tight', pad_inches=0.05,
                                dpi=400)

            else:
                warn('Matplotlib cannot be imported from Abaqus')
            np.savez(npzname, cir=cir, mer=mer, field=field)
            with open(pyname, 'w') as f:
                f.write("import os\n")
                f.write("import numpy as np\n")
                f.write("import matplotlib.pyplot as plt\n")
                f.write("tmp = np.load('abaqus_output.npz')\n")
                f.write("cir = tmp['cir']\n")
                f.write("mer = tmp['mer']\n")
                f.write("field = tmp['field']\n")
                f.write("clean = {0}\n".format(clean))
                f.write("filename = '{0}'\n".format(filename))
                f.write("plt.figure(figsize={0})\n".format(figsize))
                f.write("ax = plt.gca()\n")
                f.write("levels = np.linspace(field.min(), field.max(), {0})\n".format(
                        num_levels))
                f.write("ax.contourf(cir, mer, field, levels=levels)\n")
                f.write("ax.grid(b=None)\n")
                f.write("ax.set_aspect('{0}')\n".format(aspect))
                f.write("ax.xaxis.set_ticks_position('bottom')\n")
                f.write("ax.yaxis.set_ticks_position('left')\n")
                f.write("ax.xaxis.set_ticks([{0}, 0, {1}])\n".format(
                        -self.r2*pi, self.r2*pi))
                f.write("ax.xaxis.set_ticklabels([r'$-\pi$', '$0$', r'$+\pi$'])\n")
                f.write("ax.set_title(r'Abaqus, $PL=20 N$, $F_{{C}}=50 kN$, $w_{{PL}}=\beta$, $mm$')\n")
                f.write("if clean:\n")
                f.write("    ax.xaxis.set_ticks_position('none')\n")
                f.write("    ax.yaxis.set_ticks_position('none')\n")
                f.write("    ax.xaxis.set_ticklabels([])\n")
                f.write("    ax.yaxis.set_ticklabels([])\n")
                f.write("    ax.set_frame_on(False)\n")
                f.write("if not filename:\n")
                f.write("    filename = 'abaqus_result.png'\n")
                f.write("plt.savefig(filename, transparent=True,\n")
                f.write("            bbox_inches='tight', pad_inches=0.05, dpi=400)\n")
                f.write("plt.show()\n")
            msg('Output exported to "{0}"'.format(npzname))
            msg('Please run the file "{0}" to plot the output'.format(
                  pyname))
            return cir, mer, field
        except:
            traceback.print_exc()
            error('Opened plot could not be generated! :(')


    def SPLA(self, PLs, NLgeom=True, plot=False):
        """Runs the Single Perturbation Load Approach (SPLA)

        A set of non-linear results will be

        Parameters
        ----------
        PLs: list
            The perturbation loads used to build the knock-down curve. It must
            be a list of float values.
        NLgeom : bool, optional
            Flag passed to the ``static()`` method that tells whether a
            geometrically non-linear analysis is to be performed.

        Returns
        -------
        curves : list
            The sequence of curves, one curve for each perturbation load given
            in the input parameter ``PLs``.
            Each curve in the list is a ``dict`` object with the keys:

            =================    ==============================================
            Key                  Description
            =================    ==============================================
            ``'wall_time_s'``    The wall time for the non-linear analysis
            ``'name'``           The name of the curve. Ex: ``'PL = 1. N'``
            ``'cs'``             A list with a vector of Ritz constants for
                                 each load increment needed
            ``'increments'``     A list with the values of increments needed
            ``'wPLs'``           A list with the normal displacement at the
                                 perturbation load application point for each
                                 load increment
            ``'uTMs'``           A list containing the axial displacement for
                                 each load increment
            ``'Fcs'``            A list containing the axial reaction force
                                 for each load increment
            =================    ==============================================

        Notes
        -----
        The curves are stores in the ``ConeCyl`` parameter
        ``outputs['SPLA_curves']``.

        """
        curves = []
        for PLi, PL in enumerate(PLs):
            curve = {}
            self.forces = []
            self.add_SPL(PL)
            time1 = time.clock()
            cs = self.static(NLgeom=NLgeom)
            if plot:
                self.plot(cs[-1])
            curve['wall_time_s'] = time.clock() - time1
            curve['name'] = 'PL = {} N'.format(PL)
            curve['cs'] = cs
            curve['increments'] = self.increments
            curve['wPLs'] = []
            curve['uTMs'] = []
            curve['Fcs'] = []
            for i, c in enumerate(self.cs):
                inc = self.increments[i]
                self.uvw(c, xs=self.L/2, ts=0)
                curve['wPLs'].append(self.w[0])
                if self.pdC:
                    ts = np.linspace(0, pi*2, 1000, endpoint=False)
                    xs = np.zeros_like(ts)
                    es = self.strain(c=c, xs=xs, ts=ts, inc=inc)
                    fvec = self.F.dot(es.T)
                    Fc = -fvec[0,:].mean()*(2*self.r2*pi)
                    curve['Fcs'].append(Fc/1000)
                    curve['uTMs'].append(inc*self.uTM)
                else:
                    curve['Fcs'].append(inc*self.Fc/1000)
                    curve['uTMs'].append(c[0])
            curves.append(curve)

        self.outputs['SPLA_curves'] = curves

        return curves


    def apply_shim(self, thetadeg, width, thick, ncpts=10000):
        r"""Distributes the axial load in order to simulate a shim

        The axial load distribution `{N_{xx}}_{top}` will be adjusted such
        that the resulting displacement `u` at `x=0` (top edge) will look
        similar to a case where a shim is applied.

        Parameters
        ----------
        thetadeg : float
            Position in degrees of the center of the shim.
        width : float
            Circumferential width of the shim.
        thick : float
            Thickness of the shim.
        ncpts : int, optional
            Number of control points used in the least-squares routine.

        Returns
        -------
        ts : np.ndarray
            Positions `\theta` of the control points.
        us : np.ndarray
            Displacements `u` of the control points.

        Notes
        -----
        This function changes the ``Nxxtop`` parameter of the current
        ``ConeCyl`` object. Returning ``ts`` and ``us`` is useful for post
        processing purposes only.

        Examples
        --------
        >>> ts, us = cc.apply_shim(0., 25.4, 0.1)

        """
        ts = np.linspace(-np.pi, np.pi, ncpts)
        xs = np.zeros_like(ts)
        us = np.zeros_like(ts)
        self.static(NLgeom=False)
        thetashim = width/self.r2
        thetarad = deg2rad(thetadeg)
        theta1 = thetarad - thetashim
        theta2 = thetarad + thetashim
        uTM = self.cs[0][0]
        utop = uTM/self.cosa
        us += uTM
        shim_region = (ts >= theta1) & (ts <= theta2)
        us[shim_region] += thick
        self.fit_Nxxtop(ts, us)

        return ts, us


    def fit_Nxxtop(self, ts, us, update_Nxxtop=True):
        r"""Adjusts the axial load distribution for a desired top edge
        displacement

        Parameters
        ----------
        ts : np.ndarray
            Corrdinates `\theta` of each control point.
        us : np.ndarray
            Desired displacement `u` for each control point.
        update_Nxxtop : bool, optional
            Tells whether ``self.Nxxtop`` should be updated.

        Returns
        -------
        Nxxtop : np.ndarray
            The coefficients for the `{N_{xx}}_{top}` function.

        """
        from scipy.sparse.linalg import inv as sparseinv

        assert ts.ndim == 1 and us.ndim == 1
        assert ts.shape[0] == us.shape[0]
        xs = np.zeros_like(ts)

        if not update_Nxxtop:
            Nxxtop_backup = self.Nxxtop.copy()

        k0uuinv = sparseinv(csc_matrix(self.k0uu))

        Nxxtop_new = self.Nxxtop.copy()

        def fit_Nxxtop_residual(Nxxtop_new):
            self.Nxxtop = Nxxtop_new.copy()
            fext = self.calc_fext(silent=True)
            c = k0uuinv*fext
            self.uvw(c, xs=xs, ts=ts)
            res = (self.u - us)**2
            return res
        popt, pcov = leastsq(fit_Nxxtop_residual, x0=Nxxtop_new, maxfev=10000)

        if not update_Nxxtop:
            self.Nxxtop = Nxxtop_backup
        else:
            self.Nxxtop = popt

        return popt


    def save(self):
        """Save the ``ConeCyl`` object using ``cPickle``

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

        """
        name = self.name + '.ConeCyl'
        msg('Saving ConeCyl to {}'.format(name))
        self.analysis.calc_fext = None
        self.analysis.calc_k0 = None
        self.analysis.calc_fint = None
        self.analysis.calc_kT = None
        self._clear_matrices()

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

