DEV: development of kpanelt and analysis modules

compmech/__init__.py

@@ -1,2 +1,2 @@
-__version__ = '0.2.0'
+__version__ = '0.3.0'
 

compmech/analysis/__init__.py

@@ -0,0 +1,11 @@
+r"""
+==========================================
+Analysis Module (:mod:`compmech.analysis`)
+==========================================
+
+.. currentmodule:: compmech.analysis
+
+.. autoclass:: Analysis
+    :members:
+"""
+from analysis import Analysis

compmech/analysis/analysis.py

@@ -0,0 +1,164 @@
+import numpy as np
+from numpy import dot
+
+from compmech.sparse import solve
+from compmech.logger import msg
+from newton_raphson import _solver_NR
+from arc_length import _solver_arc_length
+
+class Analysis(object):
+    r"""Class that embodies all data required for linear/non-linear analysis
+
+    The parameters are described in the following tables:
+
+    ====================    ==================================================
+    non-linear algorithm    description
+    ====================    ==================================================
+    ``NL_method``           ``str``, ``'NR'`` for the Newton-Raphson
+                            ``'arc_length'`` for the Arc-Length method
+    ``line_search``         ``bool``, activate line_search (for Newton-Raphson
+                            methods only)
+    ``modified_NR``         ``bool``, activates the modified Newton-Raphson
+    ``compute_every_n``     ``int``, if ``modified_NR=True``, the non-linear
+                            matrices will be updated at every `n` iterations
+    ``kT_initial_state``    ``bool``, tells if the tangent stiffness matrix
+                            should be calculated already at the initial state,
+                            which is required for example when initial
+                            imperfections take place
+    ====================    ==================================================
+
+    ==============     =================================================
+    incrementation     description
+    ==============     =================================================
+    ``initialInc``     initial load increment size. In the arc-length
+                       method it will be the initial value for
+                       `\lambda`
+    ``minInc``         minimum increment size; if achieved the analysis
+                       is terminated. The arc-length method will use
+                       this parameter to terminate when the minimum
+                       arc-length increment is smaller than ``minInc``
+    ``maxInc``         maximum increment size
+    ==============     =================================================
+
+    ====================    ============================================
+    convergence criteria    description
+    ====================    ============================================
+    ``absTOL``              the convergence is achieved when the maximum
+                            residual force is smaller than this value
+    ``maxNumIter``          maximum number of iteration; if achieved the
+                            load increment is bisected
+    ``too_slow_TOL``        tolerance that tells if the convergence is too
+                            slow
+    ====================    ============================================
+
+    =====================    ===========================================
+    numerical integration    description
+    =====================    ===========================================
+    ``ni_num_cores``         number of cores used for the numerical
+                             integration
+    ``ni_method``            ``'trapz2d'`` for 2-D Trapezoidal's
+                             ``'simps2d'`` for 2-D Simpsons' integration
+    =====================    ===========================================
+
+    Parameters
+    ----------
+    calc_fext : callable, optional
+        Must return a 1-D array containing the external forces. Required for
+        linear/non-linear static analysis.
+    calc_k0 : callable, optional
+        Must return a sparse matrix containing the linear stiffness matrix.
+        Required for linear/non-linear static analysis.
+    calc_fint : callable, optional
+        Must return a 1-D array containing the internal forces. Required for
+        non-linear analysis.
+    calc_kT : callable, optional
+        Must return a sparse matrix containing the tangent stiffness matrix.
+        Required for non-linear analysis.
+
+    Returns
+    -------
+    increments : list
+        Each time increment that achieved convergence.
+    cs : list
+        The solution for each increment.
+
+    """
+
+    __slots__ = ('NL_method', 'line_search', 'modified_NR', 'compute_every_n',
+            'kT_initial_state', 'initialInc', 'minInc', 'maxInc', 'absTOL',
+            'relTOL', 'maxNumIter', 'too_slow_TOL', 'ni_num_cores',
+            'ni_method', 'nx', 'nt', 'calc_fext', 'calc_k0', 'calc_fint',
+            'calc_kT', 'increments', 'cs', 'last_analysis')
+
+
+    def __init__(self, calc_fext=None, calc_k0=None, calc_fint=None,
+            calc_kT=None):
+        # non-linear algorithm
+        self.NL_method = 'NR'
+        self.line_search = True
+        self.modified_NR = True
+        self.compute_every_n = 6
+        self.kT_initial_state = False
+        # incrementation
+        self.initialInc = 0.3
+        self.minInc = 1.e-3
+        self.maxInc = 1.
+        # convergence criteria
+        self.absTOL = 1.e-3
+        self.relTOL = 1.e-3
+        self.maxNumIter = 30
+        self.too_slow_TOL = 0.01
+        # numerical integration
+        self.ni_num_cores = 4
+        self.ni_method = 'trapz2d'
+        self.nx = 120
+        self.nt = 120
+
+        # required methods
+        self.calc_fext = calc_fext
+        self.calc_k0 = calc_k0
+        self.calc_fint = calc_fint
+        self.calc_kT = calc_kT
+
+        # outputs to be filled
+        self.increments = None
+        self.cs = None
+
+        # flag telling the last analysis
+        self.last_analysis = ''
+
+
+    def static(self, NLgeom=False, silent=False):
+        """General solver for static analyses
+
+        Selects the specific solver based on the ``NL_method`` parameter.
+
+        """
+        self.increments = []
+        self.cs = []
+
+        if NLgeom:
+            self.maxInc = max(self.initialInc, self.maxInc)
+            msg('Started Non-Linear Static Analysis', silent=silent)
+            if self.NL_method is 'NR':
+                _solver_NR(self)
+            elif self.NL_method is 'arc_length':
+                _solver_arc_length(self)
+            else:
+                raise ValueError('{0} is an invalid NL_method')
+
+        else:
+            msg('Started Linear Static Analysis', silent=silent)
+            fext = self.calc_fext()
+            k0 = self.calc_k0()
+
+            c = solve(k0, fext)
+
+            self.cs.append(c)
+            self.increments.append(1.)
+            msg('Finished Linear Static Analysis', silent=silent)
+
+        self.last_analysis = 'static'
+
+        return self.increments, self.cs
+

compmech/analysis/arc_length.py

@@ -0,0 +1,159 @@
+def _solver_arc_length(a):
+    r"""Arc-Length solver
+
+    """
+    log('Initialization...', level=1)
+    lbd_init = a.initialInc
+    lbd = lbd_init
+    last_lbd = 0.
+
+    length_inc = 1.
+    length = length_inc
+    max_length = length
+    modified_NR = a.modified_NR
+    fext = a.calc_fext(inc=1.)
+    k0 = a.calc_k0()
+    c = solve(k0, lbd*fext)
+    beta = c.dot(c)
+    fint = k0*c
+    last_fint = fint
+    step_num = 1
+    total_length = 0
+
+    if modified_NR:
+        if a.kT_initial_state:
+            log('Updating kT for initial imperfections...', level=1)
+            kT_last = a.calc_kT(c*0)
+            log('kT updated!', level=1)
+        else:
+            kT_last = k0
+        compute_NL_matrices = False
+    else:
+        compute_NL_matrices = True
+        kT_last = k0
+
+    while step_num < 100:
+        log('Attempting arc-length = {0}'.format(length), level=1)
+        min_Rmax = 1.e6
+        prev_Rmax = 1.e6
+        converged = False
+        iteration = 0
+
+        kT = kT_last
+
+        while True:
+            iteration += 1
+            log('Iteration: {}'.format(iteration), level=2)
+            if iteration > a.maxNumIter:
+                warn('Maximum number of iterations achieved!', level=2)
+                break
+
+            # applying the arc-length constraint to find the new lbd and the
+            # new c
+            dc2 = solve(kT, fext)
+            dc1 = solve(kT, (lbd*fext - fint))
+            a2 = beta + dot(dc2.T, dc2)
+            a1 = 2*beta*lbd + 2*dot(c.T, dc2) + 2*dot(dc1.T, dc2)
+            a0 = (beta*lbd**2 + dot(c.T, c) + 2*dot(c.T, dc1)
+                  + dot(dc1.T, dc1) - beta*length**2)
+            delta = a1**2 - 4*a2*a0
+            if delta<0:
+                warn('Diverged! (negative delta)', level=2)
+                break
+
+            eta = 1.
+            dlbd1 = (-a1 + np.sqrt(delta))/(2*a2)
+            dc_1 = (dc1 + dlbd1*dc2)
+            dlbd2 = (-a1 - np.sqrt(delta))/(2*a2)
+            dc_2 = (dc1 + dlbd2*dc2)
+
+            lbd = lbd + eta*dlbd1
+            c = c + eta*dc_1
+
+            if lbd > 1.:
+                warn('Diverged! (lbd > 1.)', level=2)
+                break
+
+            # computing the Non-Linear matrices
+            if compute_NL_matrices:
+                kT = a.calc_kT(c)
+            else:
+                if not modified_NR:
+                    compute_NL_matrices = True
+                #NOTE attempt to calculate fint more often than kT
+
+            fint = a.calc_fint(c)
+
+            # calculating the residual
+            Rmax = np.abs(lbd*fext - fint).max()
+            log('Rmax = {0}'.format(Rmax), level=3)
+            log('lbd = {0}'.format(lbd), level=3)
+            if Rmax <= a.absTOL:
+                converged = True
+                break
+            if (Rmax > min_Rmax and Rmax > prev_Rmax and iteration > 2):
+                warn('Diverged!', level=2)
+                break
+            else:
+                min_Rmax = min(min_Rmax, Rmax)
+            change_rate_Rmax = abs(prev_Rmax-Rmax)/abs(prev_Rmax)
+            if (iteration > 2 and change_rate_Rmax < TOO_SLOW):
+                warn('Diverged! (convergence too slow)', level=2)
+                break
+            prev_Rmax = Rmax
+
+        if converged:
+            log('Finished increment with total length = {0}'.format(
+                length), level=1)
+            a.increments.append(lbd)
+            a.cs.append(c.copy())
+            if abs(lbd - 1.) < 0.001:
+                log('Condition abs(lbd - 1.) < 0.001 achieved!', level=1)
+                break
+
+            if lbd < last_lbd and len(a.cs) > 0:
+                log('Drop of reaction load achieved!', level=1)
+                #TODO maybe when we need to detect the local snap-through
+                #     this if should be extended...
+                #break
+            last_fint = fint.copy()
+            last_lbd = lbd
+            step_num += 1
+            factor = 0.95
+            new_length = length + length_inc*factor
+            #TODO
+            #while new_length >= max_length:
+            #    new_length *= 0.9
+            log('(lambda of this step = {0})'.format(lbd), level=1)
+            log('Changing arc-length from {0} to {1}'.format(length,
+                new_length), level=1)
+            length = new_length
+            if modified_NR:
+                log('Updating kT...', level=1)
+                kT = a.calc_kT(c)
+                log('kT updated!', level=1)
+            compute_NL_matrices = False
+            kT_last = kT
+        else:
+            if len(a.cs) > 0:
+                c = a.cs[-1].copy()
+            else:
+                lbd = lbd_init
+                kT = k0
+                c = solve(kT, lbd*fext)
+            factor = 0.3 # keep in the range 0 < factor < 1.
+            fint = last_fint
+            lbd = last_lbd
+            max_length = max(length, max_length)
+            old_length = length
+            length -= length_inc
+            length_inc *= factor
+            length += length_inc
+            log('Diverged - reducing arc-length from {0} to {1}'.format(
+                old_length, length), level=1)
+            if length_inc < a.minInc:
+                log('Minimum arc-length achieved!', level=1)
+                break
+
+    log('Finished Non-Linear Static Analysis')
+    log('with a total arc-length {0}'.format(length), level=1)