Commits

Panagiotis Mavrogiorgos committed 840455f Merge

merged from default

Comments (0)

Files changed (16)

 syntax: glob
 
+*.o
+*.c
 *.pyc
-*.wpr
+*.so
+*.wpr
+*.txt
+*.tex
+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+# module : constants.py
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+"""
+constants.py Module
+
+This module contains the constants of the program.
+"""
+
+##------------------------------------------------------------------------------
+# Importing necessary libraries
+##------------------------------------------------------------------------------
+from __future__ import division
+
+##------------------------------------------------------------------------------
+# Defining constants
+##------------------------------------------------------------------------------
+## The maximum number of iterations allowed in order to achieve convergence
+MAX_ITERATIONS = 500
+## The tolerance for which convergence is assumed.
+TOLERANCE = 1E-20
+## The length of the infinitestimal element in meters. 1/100 is one cm.
+## For development 1/100 is suggested. For final analysis 1/1000.
+INFINITESIMAL_LENGTH = 1/10
+## the length (in m) where we suppose an omoiothetic decrease of the bond-slip law.
+LENGTH_OF_OMOIOTHECY = 0.05
+## If True then omoiothecy length is taken under consideration for the front of the specimen.
+OMOIOTHECY_FRONT = True
+#OMOIOTHECY_FRONT = False
+## If True then omoiothecy length is taken under consideration for the back of the specimen.
+#OMOIOTHECY_BACK = True
+OMOIOTHECY_BACK = False
+## There are numerical instabilities when unloading to zero, so we stop at +-
+## LOAD_STABILITY.
+LOAD_STABILITY = 0.1
+##------------------------------------------------------------------------------
 #! /usr/bin/env python2
 # -*- coding: utf-8 -*-
-# module : new.py
+# module : bond.py
 ##------------------------------------------------------------------------------
 
 ##------------------------------------------------------------------------------
 # Licence : GPL v3
 ##------------------------------------------------------------------------------
 
+"""
+Contains the classes with the Bond definition
+"""
+
 from __future__ import division
 from __future__ import print_function
-from copy import deepcopy as dcopy
+from __future__ import unicode_literals
+
+from copy import copy as ccopy
 from math import sqrt, exp
-from matplotlib.pyplot import (figure, axhline, axvline, title, legend,
-                               show, xlabel, ylabel, grid, plot,
-                               close as closefig)
 import numpy
 from numpy import (float64 as np_float,
                    isnan,
                    where)
-from concrete import ModelCode2010
 from lines import (line_slope, line_x, line_y, isascending, isdescending,
                    interpolation, calc_area, line_intersection)
-from bunch import Bunch
 
 class Bond(object):
     """ Abstract class"""
 
         return t
 
-    def db_initialize_bond(*args, **kwargs):
+    def initialize_db(*args, **kwargs):
         """ Must be overriden in the derived classes"""
         raise NotImplementedError
 
-    def db_store_bond(*args, **kwargs):
+    def store_to_db(*args, **kwargs):
         """ Must be overriden in the derived classes"""
         raise NotImplementedError
 
         tmax = fck / 3
         tr = friction_factor * tmax
 
+        #slip_mono = numpy.array([0, 0.005,  0.05,  0.2,   1,   2, 10, 15])
+        #bond_mono = numpy.array([0, 3.133, 6.267,  9.4, 9.4, 9.4,  5,  5])
+
         slip_mono = [0, 0.01,     0.1,  0.5,    1, 10.5, 15]
         bond_mono = [0,  fct, 2 * fct, tmax, tmax,   tr, tr]
 
         # calculate slopes of the monotonic curve segments.
         self.slope1 = line_slope(slip_mono[0], bond_mono[0],
                                  slip_mono[1], bond_mono[1])
-        self.slope2 = line_slope(slip_mono[1], bond_mono[1],
-                                 slip_mono[2], bond_mono[2])
-        self.slope3 = line_slope(slip_mono[2], bond_mono[2],
-                                 slip_mono[3], bond_mono[3])
 
     def _calc_reduced_monotonic_bond(self, n):
         """
         Only bond values are retured. No sign is assigned.
         """
         if n <= 10:
-            bond = (1 - sqrt((n//2) / 10)) * self.bond_mono
+            bond = (1 - sqrt((1 + n//2) / 10)) * self.bond_mono
         else:
             bond = (1 - sqrt(5 / 10)) * self.bond_mono
         return bond
         return tf
 
     def _bond_unload(self, s, t, smin, smax, n):
+        """
+        Calculates the unloading slip-bond constitutive law.
+        """
         # Check for errors
         assert smin <= 0, "smin can't be positive!!!"
         assert smax >= 0, "smax can't be negative!!!"
 
+        slope1 = self.slope1
+
         # get reduced monotonic curves and assign sign
         red_bond = -self._calc_reduced_monotonic_bond(n)
         red_slip = -self.slip_mono
         # tmin is the bond that corresponds to smin and lies on red_bond.
         tf = self._calc_tf(t, n)
 
-        # If during the reloading, slip sign hasn't changed (aka it is negative)
-        # then depending on the bond value the first points of the bond-slip
-        # curve are calculated.
+        # If during the reloading, slip sign hasn't changed (eg it is negative)
+        # then, the first points of the bond-slip
+        # curve are calculated. Depending on the bond value there are 3 cases.
+        # a) bond lower than zero
+        # b) bond equal to zero
+        # c) bond greter than zero.
+        #
+        # a) The 1st point is the intersection of the x axis and the lines that
+        # passes through (s, t) with slope equal to slope1 (e.g. (s', 0)).
+        # The 2nd point is (s, t) . If there has already been slipping in the
+        # negative direction, the 3rd point is (smin, tmin).
+        #
+        # b) The 1st point is (s, t) (eg (s, 0)). If there has been slipping
+        # in the negative direaction, the 2nd point is (smin, tmin).
+        #
+        # c) The 1st point is (s, t). The second point is the intersection of
+        # the x-axis and the line that passes through (s, t) with slope equal
+        # to slope1 (e.g. (s', 0)). If there has been slipping in the negative
+        # direction in the previous cycles, then the 3rd point is (smin, tmin)
         if s <= 0:
             if t < 0:
-                slip = [line_x(self.slope1, s, t, 0), s]
+                slip = [line_x(slope1, s, t, 0), s]
                 bond = [0, t]
                 if smin < 0 and smin < s:
                     slip.append(smin)
                     slip.append(smin)
                     bond.append(tmin)
             else:
-                slip = [s, line_x(self.slope1, s, t, 0)]
+                slip = [s, line_x(slope1, s, t, 0)]
                 bond = [t, 0]
                 if smin < 0 and smin < s:
                     slip.append(smin)
                     bond.append(tmin)
-        # If slipping has occured but it is very small (smax < 0.01mm),
-        # then elastic behavior is assumed untl slip changes sign.
-        # Then linear behavior up to (smin, tmin).
+
+        # If during the reloading, slip sign has changed (eg it became positive)
+        # but the slip value is very small, (s < 0.01mm), then it is assumed
+        # that an elastic behavior is followed, until the next change of sign.
         elif s < self.slip_mono[1]:
             slip = [s, 0]
             bond = [t, 0]
-        # If slipping is greater than smax<0.01 but small enough that unloading
-        # with slope1 results in a change of slip sign, then the following rule
-        # is followed
-        elif line_x(self.slope1, s, t, tf) <= 0:
-            tf = line_y(self.slope1, s, t, 0)
+
+        # If slipping is greater than (s > 0.01mm) but small enough that
+        # unloading with slope1 results in a change of slip sign before the bond
+        # value reaches tf, then the 1st point is (s, t).
+        # The 2nd point is the intersection of the x_axis and the line that
+        # passes through (s, t) with slope equal to slope1 (eg (s', 0)).
+        # The 3rd point is the intersection of the y-axis and the line that
+        # passes through (s, t) with slope equal to slope1 (eg (0, t')).
+        elif line_x(slope1, s, t, tf) <= 0:
+            tf = line_y(slope1, s, t, 0)
             slip = [s,
-                    line_x(self.slope1, s, t, 0.),
+                    line_x(slope1, s, t, 0.),
                     0.]
             bond = [t, 0, tf]
-        # If slipping is significant then normal unloading is followed.
+        # If slipping is significant then the first point is (s, t)
+        # The 2nd point is the intersection of the x-axis with the line that
+        # passes through (s, t) with slope equal to slope1 (eg (s', 0)).
+        # The 3rd point is the point that belongs to the line that passes through
+        # (s, t) with slope equal to slope1 and it's ordinate is equal to tf
+        # The 4th point is (0, tf)
+        #
+        # In case there has been significant slipping in the opposite direction,
+        # (eg smin < 0.01mm) then the 5th point is (smin/2, tf).
+        # The 6th point's abscissa is the point with abscissa equal to the min
+        # value of smin and the abscissa of the reduced monotonic curve with
+        # bond value equal to tf.
+        # The 6th's point ordinate is the ordinate of the reduced monotonic
+        # curve's point with slip value equal to the 6th point's abscissa.
         else:
             slip = [s,
-                    line_x(self.slope1, s, t, 0.),
-                    line_x(self.slope1, s, t, tf),
+                    line_x(slope1, s, t, 0.),
+                    line_x(slope1, s, t, tf),
                     0.]
             bond = [t, 0, tf, tf]
-            # If no slipping in the opposite direction has occured
-            # then the reduced monotonic curve is followed.
-            # If slipping has occured then the unloading rules are followed.
             if smin < red_slip[1]:
                 slip.append(smin / 2)
                 bond.append(tf)
                 slip.append(min(smin, interpolation(red_bond, red_slip, tf)))
                 bond.append(interpolation(red_slip, red_bond, slip[-1]))
 
-        # Append the values of the reduced envelope
+        #After the calculation of the first points, the rest of the points belong
+        # to the reduced envelope.
         for j, value in enumerate(red_slip):
             if value < slip[-1]:
                 slip.append(red_slip[j])
                 bond.append(red_bond[j])
 
         # Sanity check
-        # If values are ok, return them, else raise an error.
         assert isdescending(slip), "Slip array must be monotonous"
         assert len(slip) == len(bond), "Slip - bond arrays must have the same size!!"
         assert not (any(isnan(slip)) or any(isnan(bond))), "Nan values!!!"
 
         return slip, bond
 
-    def define_bond_unloading(self, slip, bond, s, t, smin, smax, N):
+    def define_bond_unloading(self, s, t, smin, smax, N):
         """
         Defines bond's unloading constitutive laws.
         Should be called everytime there is a change in loading sign.
         """
+        slip, bond = self._bond_unload(s, t, smin, smax, N)
 
-        smax = max(smax, s)
+        return slip, bond
 
-        slip_new, bond_new = self._bond_unload(s, t, smin, smax, N)
-
-        return slip_new, bond_new, smax
-
-    def define_bond_reloading(self, slip, bond, s, t, smin, smax, N):
+    def define_bond_reloading(self, s, t, smin, smax, N):
         """
         Defines bond's reloading constitutive laws.
         Should be called everytime there is a change in loading sign.
         """
-        smin = min(smin, s)
+        slip, bond = self._bond_reload(s, t, smin, smax, N)
 
-        slip_new, bond_new = self._bond_reload(s, t, smin, smax, N)
-
-        return slip_new, bond_new, smin
+        return slip, bond
 
     @staticmethod
-    def db_initialize_bond(db):
-        assert type(db) is Bunch, "db must be a bunch instance!!!"
+    def initialize_db(db):
 
-        db.slip_tot = []
-        db.bond_tot = []
-        db.smin_tot = []
-        db.tmin_tot = []
-        db.smax_tot = []
-        db.tmax_tot = []
+        db["slip_tot"] = []
+        db["bond_tot"] = []
+        db["smin_tot"] = []
+        db["tmin_tot"] = []
+        db["smax_tot"] = []
+        db["tmax_tot"] = []
 
     @staticmethod
-    def db_store_bond(db, slip, bond, smin, smax):
-        assert type(db) is Bunch, "db must be a bunch instance!!!"
+    def store_to_db(db, slip, bond, smin, smax):
 
-        db.slip_tot.append(dcopy(slip))
-        db.bond_tot.append(dcopy(bond))
-        db.smin_tot.append(dcopy(smin))
-        #db.tmin_tot.append(dcopy(tmin))
-        db.smax_tot.append(dcopy(smax))
-        #db.tmax_tot.append(dcopy(tmax))
-
-class TestTassios(Tassios):
-    def __init__(self, fck, fctm):
-        super(TestTassios, self).__init__(fck, fctm)
-
-        self.db = Bunch()
-
-    def plot_bond(self, numbers, x_label="Slip (mm)", y_label = "Bond Stress (MPa)"):
-        """
-        Plots the local bond - local slip history for node i.
-        """
-        figure(0)
-        title("Tassios's bond-slip law")
-        plot(self.slip_mono, self.bond_mono, label="Pos Mono")
-        plot(-self.slip_mono, -self.bond_mono, label="Neg Mono")
-
-        for j in numbers:
-            plot(self.db.slip_tot[j], self.db.bond_tot[j], label=str(j))
-        axhline(y=0, color = "black")
-        axvline(x=0, color = "black")
-        ylabel(y_label)
-        xlabel(x_label)
-        legend(loc = "best")
-        grid(True)
-
-        show()
-
-    def recycling_tests_positive(self, s, times):
-        """
-        s : float
-        times : int
-        """
-        db = self.db
-        self.db_initialize_bond(db)
-
-        # Monotonic pull-out
-        N = 0
-        smin = 0; tmin = 0; smax = 0; tmax = 0;
-        slip = self.slip_mono
-        bond = self.bond_mono
-
-        self.db_store_bond(db, slip, bond, smin, smax)
-
-        # Start the cycles
-        for i in range(0, times):
-            N += 1
-            t = self.get_bond(slip, bond, s)
-            slip, bond, smax = self.define_bond_unloading(
-                slip, bond, s, t, smin, smax, N)
-            self.db_store_bond(db, slip, bond, smin, smax)
-
-            N += 1
-            t = self.get_bond(slip, bond, -s)
-            slip, bond, smin = self.define_bond_reloading(
-                slip, bond, -s, t, smin, smax, N)
-            self.db_store_bond(db, slip, bond, smin, smax)
-
-        # End with a pull-out
-        N += 1
-        t = self.get_bond(slip, bond, s)
-        slip, bond, smax = self.define_bond_unloading(
-            slip, bond, s, t, smin, smax, N)
-        self.db_store_bond(db, slip, bond, smin, smax)
-
-
-    def recycling_tests_negative(self, s, times):
-        """
-        s : float
-        times : int
-        """
-        db = self.db
-        self.db_initialize_bond(db)
-
-        # Monotonic pull-out
-        N = 0
-        smin = 0; tmin = 0; smax = 0; tmax = 0;
-        slip = -self.slip_mono
-        bond = -self.bond_mono
-
-        self.db_store_bond(self.db, slip, bond, smin, smax)
-
-        # Start the cycles
-        for i in range(0, times):
-            N += 1
-            t = self.get_bond(slip, bond, -s)
-            slip, bond, smin = self.define_bond_reloading(
-                slip, bond, -s, t, smin, smax, N)
-            self.db_store_bond(self.db, slip, bond, smin, smax)
-
-            N += 1
-            t = self.get_bond(slip, bond, s)
-            slip, bond, smax = self.define_bond_unloading(
-                slip, bond, s, t, smin, smax, N)
-            self.db_store_bond(self.db, slip, bond, smin, smax)
-
-        # End with a push in
-        N += 1
-        t = self.get_bond(slip, bond, -s)
-        slip, bond, smin = self.define_bond_reloading(
-            slip, bond, -s, t, smin, smax, N)
-        self.db_store_bond(self.db, slip, bond, smin, smax)
-
-    def repeating_tests_positive(self, s, times):
-        """
-        Repeating loading, start from positive slips
-
-        s : float
-        times : int
-        """
-        db = self.db
-        self.db_initialize_bond(db)
-
-        # Monotonic pull-out
-        N = 0
-        smin = 0; tmin = 0; smax = 0; tmax = 0;
-        slip = self.slip_mono
-        bond = self.bond_mono
-
-        self.db_store_bond(db, slip, bond, smin, smax)
-
-        # Start the cycles
-        for i in range(0, times):
-            N += 1
-            t = self.get_bond(slip, bond, s)
-            slip, bond, smax = self.define_bond_unloading(
-                slip, bond, s, t, smin, smax, N)
-            self.db_store_bond(db, slip, bond, smin, smax)
-
-            N += 1
-            t = self.get_bond(slip, bond, 0)
-            slip, bond, smin = self.define_bond_reloading(
-                slip, bond, 0, t, smin, smax, N)
-            self.db_store_bond(db, slip, bond, smin, smax)
-
-        # End with a pull-out
-        N += 1
-        t = self.get_bond(slip, bond, s)
-        slip, bond, smax = self.define_bond_unloading(
-            slip, bond, s, t, smin, smax, N)
-        self.db_store_bond(db, slip, bond, smin, smax)
-
-    def repeating_tests_negative(self, s, times):
-        """
-        Repeating loading, start from negative slips
-
-        s : float
-        times : int
-        """
-        db = self.db
-        self.db_initialize_bond(db)
-
-        # Monotonic pull-out
-        N = 0
-        smin = 0; tmin = 0; smax = 0; tmax = 0;
-        slip = -self.slip_mono
-        bond = -self.bond_mono
-
-        self.db_store_bond(db, slip, bond, smin, smax)
-
-        # Start the cycles
-        for i in range(0, times):
-            N += 1
-            t = self.get_bond(slip, bond, -s)
-            slip, bond, smin = self.define_bond_reloading(
-                slip, bond, -s, t, smin, smax, N)
-            self.db_store_bond(db, slip, bond, smin, smax)
-
-            N += 1
-            t = self.get_bond(slip, bond, 0)
-            slip, bond, smax = self.define_bond_unloading(
-                slip, bond, 0, t, smin, smax, N)
-            self.db_store_bond(db, slip, bond, smin, smax)
-
-        # End with a pull-out
-        N += 1
-        t = self.get_bond(slip, bond, -s)
-        slip, bond, smin = self.define_bond_reloading(
-            slip, bond, -s, t, smin, smax, N)
-        self.db_store_bond(db, slip, bond, smin, smax)
+        db["slip_tot"].append(ccopy(slip))
+        db["bond_tot"].append(ccopy(bond))
+        db["smin_tot"].append(ccopy(smin))
+        #db["tmin_tot"].append(ccopy(tmin))
+        db["smax_tot"].append(ccopy(smax))
+        #db["tmax_tot"].append(ccopy(tmax))
 
 class Eligenhausen(Bond):
     def __init__(self, fck, fctm):
         slip[-1] = slip[-2]
         bond[-1] = 0
 
-        E0 = calc_area(slip, bond)
+        # Converting slip and bond arrays to lists because
+        # calc_area may be cythonized and it will require lists
+        # more info at lines.pxd
+        E0 = calc_area(slip.tolist(), bond.tolist())
         E0f = slip[-2] * bond[-2]
 
         self.E0 = abs(E0)
         # calculate slopes of the monotonic curve segments.
         self.slope1 = line_slope(slip_mono[0], bond_mono[0],
                                  slip_mono[1], bond_mono[1])
-        self.slope2 = line_slope(slip_mono[1], bond_mono[1],
-                                 slip_mono[2], bond_mono[2])
-        self.slope3 = line_slope(slip_mono[2], bond_mono[2],
-                                 slip_mono[3], bond_mono[3])
 
     def _bond_unload(self, s, t, smin, smax, tmin, tf, d, g):
         """
         # get reduced monotonic curves and assign sign
         red_bond = -self._calc_reduced_monotonic_bond(d, g)
         red_slip = -self.slip_mono
+        tmin = interpolation(red_slip, red_bond, smin)
 
         # If during the reloading, slip sign hasn't changed (aka it is negative)
         # then depending on the bond value the first points of the bond-slip
         # curve are calculated.
         if s <= 0:
-            tmin = interpolation(red_slip, red_bond, smin)
             if t < 0:
-                slip = [line_x(self.slope1, s, t, 0), s, smin]
-                bond = [0, t, tmin]
+                slip = [line_x(self.slope1, s, t, 0), s]
+                bond = [0, t]
+                if smin < 0 and smin < s:
+                    slip.append(smin)
+                    bond.append(tmin)
             elif t == 0:
-                slip = [s, smin]
-                bond = [t, tmin]
+                slip = [s]
+                bond = [t]
+                if smin < 0 and smin < s:
+                    slip.append(smin)
+                    bond.append(tmin)
             else:
-                slip = [s, line_x(self.slope1, s, t, 0), smin]
-                bond = [t, 0, tmin]
+                slip = [s, line_x(self.slope1, s, t, 0)]
+                bond = [t, 0]
+                if smin < 0 and smin < s:
+                    slip.append(smin)
+                    bond.append(tmin)
         # If slipping has occured but it is very small (smax < 0.01mm),
         # then elastic behavior is assumed.
         elif s < self.slip_mono[1]:
             #
             # Calculate the point of intersection of the unloading segment of
             # the previous half-circle and the frictional part.
-            temp_x, temp_y = line_intersection(0, slip[-1], tf, self.slope1, smin, tmin)
-            assert temp_y == tf # just for checking.
-            slip.append(temp_x)
-            bond.append(tf)
+            temp_x, temp_y = line_intersection(0, slip[-1], tf,
+                                               self.slope1, smin, tmin)
+            if temp_x < slip[-1]:
+                slip.append(temp_x)
+                bond.append(tf)
             # I don't remember why slip is calculated like that...
             # Propably it is for cases where smax is very small.
             # It shouldn't really matter. Propably, just smooths. the curve.
             # Now that (smin/2) is no longer used, it should really matter.
-            slip.append(min(smin, interpolation(red_bond, red_slip, tf)))
-            bond.append(interpolation(red_slip, red_bond, slip[-1]))
+            temp = min(smin, interpolation(red_bond, red_slip, tf))
+            if temp < temp_x:
+                slip.append(temp)
+                bond.append(interpolation(red_slip, red_bond, temp))
+            #slip.append(min(smin, interpolation(red_bond, red_slip, tf)))
+            #bond.append(interpolation(red_slip, red_bond, slip[-1]))
 
         ## If no slipping in the opposite direction has occured
         ## then the reduced monotonic curve is followed.
         # get reduced monotonic curves and assign sign
         red_bond = self._calc_reduced_monotonic_bond(d, g)
         red_slip = self.slip_mono
+        tmax = interpolation(red_slip, red_bond, smax)
 
         # If during the unloading, slip sign hasn't changed (aka it is positive)
         # then depending on the bond value the first points of the bond-slip
         # curve are calculated.
         if s >= 0:
-            tmax = interpolation(red_slip, red_bond, smax)
             if t > 0:
-                slip = [line_x(self.slope1, s, t, 0), s, smax]
-                bond = [0, t, tmax]
+                slip = [line_x(self.slope1, s, t, 0), s]
+                bond = [0, t]
+                if smax > 0 and smax > s:
+                    slip.append(smax)
+                    bond.append(tmax)
             elif t == 0:
-                slip = [s, smax]
-                bond = [t, tmax]
+                slip = [s]
+                bond = [t]
+                if smax > 0 and smax > s:
+                    slip.append(smax)
+                    bond.append(tmax)
             else:
-                slip = [s, line_x(self.slope1, s, t, 0), smax]
-                bond = [t, 0, tmax]
+                slip = [s, line_x(self.slope1, s, t, 0)]
+                bond = [t, 0]
+                if smax > 0 and smax > s:
+                    slip.append(smax)
+                    bond.append(tmax)
+
         # If slipping has occured but it is very small (smin < 0.01)
         # then elastic behavior is assumed.
         elif s > -self.slip_mono[1]:
             #
             # Calculate the point of intersection of the unloading segment of
             # the previous half-circle and the frictional part.
-            temp_x, temp_y = line_intersection(0, slip[-1], tf, self.slope1, smax, tmax)
-            assert temp_y == tf # just for checking.
-            slip.append(temp_x)
-            bond.append(tf)
+            temp_x, temp_y = line_intersection(0, slip[-1], tf,
+                                               self.slope1, smax, tmax)
+            if temp_x > slip[-1]:
+                slip.append(temp_x)
+                bond.append(tf)
             # I don't remember why slip is calculated like that...
             # Propably it is for cases where smax is very small.
             # It shouldn't really matter. Propably, just smooths. the curve.
             # Now that (smax/2) is no longer used, it should really matter.
-            slip.append(max(smax, interpolation(red_bond, red_slip, tf)))
-            bond.append(interpolation(red_slip, red_bond, slip[-1]))
+            temp = max(smax,interpolation(red_bond, red_slip, tf))
+            if temp > temp_x:
+                slip.append(temp)
+                bond.append(interpolation(red_slip, red_bond, temp))
+            #temp_x, temp_y = line_intersection(0, slip[-1], tf, self.slope1, smax, tmax)
+            #slip.append(max(smax, temp_x))
+            #bond.append(interpolation(red_slip, red_bond, slip[-1]))
         # Append the values of the reduced envelope
         for j, value in enumerate(red_slip):
             if value > slip[-1]:
         Ef = self._calc_Ef(s, smin, smax, tf)
         Efsum += Ef
 
-        slip_new, bond_new = self._bond_unload(s, t, smin, smax, tmin, tf, d, g)
-        self.db_store_bond(self.db, slip, bond, smin, tmin, smax, tmax,
-                           E1, E2, Esum, d, g, df, tf, Ef, Efsum)
-        return slip_new, bond_new, smin, smax, tmax, Esum, Efsum
+        slip, bond= self._bond_unload(s, t, smin, smax, tmin, tf, d, g)
+
+        return smax, tmax, E1, E2, Esum, d, g, df, tf, Ef, Efsum, slip, bond
 
     def define_bond_reloading(self, s, t, slip, bond,
                               smin, smax, tmax, Esum, Efsum):
         Ef = self._calc_Ef(s, smin, smax, tf)
         Efsum += Ef
 
-        slip_new, bond_new = self._bond_reload(s, t, smin, smax, tmax, tf, d, g)
-        self.db_store_bond(self.db, slip, bond, smin, tmin, smax, tmax,
-                           E1, E2, Esum, d, g, df, tf, Ef, Efsum)
-        return slip_new, bond_new, smin, smax, tmin, Esum, Efsum
+        slip, bond = self._bond_reload(s, t, smin, smax, tmax, tf, d, g)
+
+        return smin, tmin, E1, E2, Esum, d, g, df, tf, Ef, Efsum, slip, bond
 
     @staticmethod
-    def db_initialize_bond(db):
-        assert type(db) is Bunch, "db must be a Bunch instance!!!"
+    def initialize_db(db):
 
-        db.smin_tot = []
-        db.tmin_tot = []
-        db.smax_tot = []
-        db.tmax_tot = []
-        db.slip_tot = []
-        db.bond_tot = []
+        db["smin_tot"]= []
+        db["tmin_tot"] = []
+        db["smax_tot"] = []
+        db["tmax_tot"] = []
+        db["slip_tot"] = []
+        db["bond_tot"] = []
 
-        db.E1_tot = []
-        db.E2_tot = []
-        db.Esum_tot = []
-        db.d_tot = []
-        db.g_tot = []
-        db.Ef_tot = []
-        db.Efsum_tot = []
-        db.df_tot = []
-        db.tf_tot = []
+        db["E1_tot"] = []
+        db["E2_tot"] = []
+        db["Esum_tot"] = []
+        db["d_tot"] = []
+        db["g_tot"] = []
+        db["Ef_tot"] = []
+        db["Efsum_tot"] = []
+        db["df_tot"] = []
+        db["tf_tot"] = []
 
     @staticmethod
-    def db_store_bond(db, slip, bond, smin, tmin, smax, tmax,
+    def store_to_db(db, slip, bond, smin, tmin, smax, tmax,
                       E1, E2, Esum, d, g, Ef, Efsum, df, tf):
 
-        assert type(db) is Bunch, "db must be a bunch instance!!!"
+        db["slip_tot"].append(ccopy(slip))
+        db["bond_tot"].append(ccopy(bond))
+        db["smin_tot"].append(ccopy(smin))
+        db["tmin_tot"].append(ccopy(tmin))
+        db["smax_tot"].append(ccopy(smax))
+        db["tmax_tot"].append(ccopy(tmax))
 
-        db.slip_tot.append(dcopy(slip))
-        db.bond_tot.append(dcopy(bond))
-        db.smin_tot.append(dcopy(smin))
-        db.tmin_tot.append(dcopy(tmin))
-        db.smax_tot.append(dcopy(smax))
-        db.tmax_tot.append(dcopy(tmax))
+        db["E1_tot"].append(ccopy(E1))
+        db["E2_tot"].append(ccopy(E2))
+        db["Esum_tot"].append(ccopy(Esum))
+        db["d_tot"].append(ccopy(d))
+        db["g_tot"].append(ccopy(g))
+        db["Ef_tot"].append(ccopy(Ef))
+        db["Efsum_tot"].append(ccopy(Efsum))
+        db["df_tot"].append(ccopy(df))
+        db["tf_tot"].append(ccopy(tf))
 
-        db.E1_tot.append(dcopy(E1))
-        db.E2_tot.append(dcopy(E2))
-        db.Esum_tot.append(dcopy(Esum))
-        db.d_tot.append(dcopy(d))
-        db.g_tot.append(dcopy(g))
-        db.Ef_tot.append(dcopy(Ef))
-        db.Efsum_tot.append(dcopy(Efsum))
-        db.df_tot.append(dcopy(df))
-        db.tf_tot.append(dcopy(tf))
-
-class TestEligenhausen(Eligenhausen):
-    """Tests Eligenhausen's bond-slip law"""
-    def __init__(self, fck, fctm):
-        """"""
-        super(TestEligenhausen, self).__init__(fck, fctm)
-
-        self.db = Bunch()
-        self.db_initialize_bond(self.db)
-
-    def plot_bond(self, numbers, x_label="Slip (mm)", y_label = "Bond Stress (MPa)"):
-        """
-        Plots the local bond - local slip history for node i.
-        """
-        figure(0)
-        title("Eligenhausen's bond-slip law")
-        plot(self.slip_mono, self.bond_mono, label="Pos Mono")
-        plot(-self.slip_mono, -self.bond_mono, label="Neg Mono")
-
-        for j in numbers:
-            plot(self.db.slip_tot[j], self.db.bond_tot[j], label=str(j))
-        axhline(y=0, color = "black")
-        axvline(x=0, color = "black")
-        ylabel(y_label)
-        xlabel(x_label)
-        legend(loc = "best")
-        grid(True)
-
-        show()
-
-    def run_tests(self, s, times):
-        """
-        s : float
-        times : int
-        """
-        # Monotonic pull-out
-        N = 0
-        smin = 0; tmin = 0; smax = 0; tmax = 0;
-        E1 = 0; E2 = 0; Ef = 0; Esum = 0; Efsum = 0
-        d = 0; g = 0; df = 0; tf = 0
-        slip = self.slip_mono
-        bond = self.bond_mono
-
-        self.db_store_bond(self.db, slip, bond, smin, tmin, smax, tmax,
-                           E1, E2, Esum, d, g, df, tf, Ef, Efsum)
-
-        # Start the cycles
-        for i in range(0, times):
-            N += 1
-            t = self.get_bond(slip, bond, s)
-            slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_unloading(
-                s, t, slip, bond, N, smin, smax, tmin, Esum, Efsum)
-
-            N += 1
-            t = self.get_bond(slip, bond, -s)
-            slip, bond, smin, smax, tmin, Esum, Efsum = self.define_bond_reloading(
-                -s, t, slip, bond, N, smin, smax, tmax, Esum, Efsum)
-
-        # end with a pull-out
-        N += 1
-        t = self.get_bond(slip, bond, s)
-        slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_unloading(
-            s, t, slip, bond, N, smin, smax, tmin, Esum, Efsum)
+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : bond_test.py
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+"""
+Defines classes for testing the Tassios and Eligenausen's slip-bond models
+"""
+
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+from matplotlib.pyplot import (figure, axhline, axvline, title, legend,
+                               show, xlabel, ylabel, grid, plot,
+                               close as closefig)
+
+from concrete import ModelCode2010
+from bond import Tassios
+from bond import Eligenhausen
+
+class TestTassios(Tassios):
+    def __init__(self, fck, fctm):
+        super(TestTassios, self).__init__(fck, fctm)
+
+        self.db = {}
+
+    def plot_bond(self, numbers, x_label="Slip (mm)", y_label = "Bond Stress (MPa)"):
+        """
+        Plots the local bond - local slip history for node i.
+        """
+        figure(0)
+        title("Tassios's bond-slip law")
+        plot(self.slip_mono, self.bond_mono, label="Pos Mono")
+        plot(-self.slip_mono, -self.bond_mono, label="Neg Mono")
+
+        for j in numbers:
+            plot(self.db["slip_tot"][j], self.db["bond_tot"][j], label=str(j))
+        axhline(y=0, color = "black")
+        axvline(x=0, color = "black")
+        ylabel(y_label)
+        xlabel(x_label)
+        legend(loc = "best")
+        grid(True)
+
+        show()
+
+    def recycling_tests_positive(self, s, times):
+        """
+        s : float
+        times : int
+        """
+        db = self.db
+        self.initialize_db(db)
+
+        # Monotonic pull-out
+        N = 0
+        smin = 0; smax = 0
+        slip = self.slip_mono
+        bond = self.bond_mono
+
+        self.store_to_db(db, slip, bond, smin, smax)
+
+        # Start the cycles
+        for i in range(0, times):
+            # Pull-out
+            N += 1
+            t = self.get_bond(slip, bond, s)
+            smax = max(smax, s)
+
+            slip, bond = self.define_bond_unloading(
+                s, t, smin, smax, N)
+            self.store_to_db(db, slip, bond, smin, smax)
+
+            # Push-in
+            N += 1
+            t = self.get_bond(slip, bond, -s)
+            smin = min(smin, s)
+
+            slip, bond = self.define_bond_reloading(
+                -s, t, smin, smax, N)
+
+            self.store_to_db(db, slip, bond, smin, smax)
+
+        # End with a pull-out
+        N += 1
+        t = self.get_bond(slip, bond, s)
+        smax = max(smax, s)
+
+        slip, bond = self.define_bond_unloading(
+            s, t, smin, smax, N)
+
+        self.store_to_db(db, slip, bond, smin, smax)
+
+
+    def recycling_tests_negative(self, s, times):
+        """
+        s : float
+        times : int
+        """
+        db = self.db
+        self.initialize_db(db)
+
+        # Monotonic pull-out
+        N = 0
+        smin = 0; smax = 0
+        slip = -self.slip_mono
+        bond = -self.bond_mono
+
+        self.store_to_db(self.db, slip, bond, smin, smax)
+
+        # Start the cycles
+        for i in range(0, times):
+            N += 1
+            t = self.get_bond(slip, bond, -s)
+            slip, bond, smin = self.define_bond_reloading(
+                slip, bond, -s, t, smin, smax, N)
+            self.store_to_db(self.db, slip, bond, smin, smax)
+
+            N += 1
+            t = self.get_bond(slip, bond, s)
+            slip, bond, smax = self.define_bond_unloading(
+                slip, bond, s, t, smin, smax, N)
+            self.store_to_db(self.db, slip, bond, smin, smax)
+
+        # End with a push in
+        N += 1
+        t = self.get_bond(slip, bond, -s)
+        slip, bond, smin = self.define_bond_reloading(
+            slip, bond, -s, t, smin, smax, N)
+        self.store_to_db(self.db, slip, bond, smin, smax)
+
+    def repeating_tests_positive(self, s, times):
+        """
+        Repeating loading, start from positive slips
+
+        s : float
+        times : int
+        """
+        db = self.db
+        self.initialize_db(db)
+
+        # Monotonic pull-out
+        N = 0
+        smin = 0; smax = 0
+        slip = self.slip_mono
+        bond = self.bond_mono
+
+        self.store_to_db(db, slip, bond, smin, smax)
+
+        # Start the cycles
+        for i in range(0, times):
+            N += 1
+            t = self.get_bond(slip, bond, s)
+            slip, bond, smax = self.define_bond_unloading(
+                slip, bond, s, t, smin, smax, N)
+            self.store_to_db(db, slip, bond, smin, smax)
+
+            N += 1
+            t = self.get_bond(slip, bond, 0)
+            slip, bond, smin = self.define_bond_reloading(
+                slip, bond, 0, t, smin, smax, N)
+            self.store_to_db(db, slip, bond, smin, smax)
+
+        # End with a pull-out
+        N += 1
+        t = self.get_bond(slip, bond, s)
+        slip, bond, smax = self.define_bond_unloading(
+            slip, bond, s, t, smin, smax, N)
+        self.store_to_db(db, slip, bond, smin, smax)
+
+    def repeating_tests_negative(self, s, times):
+        """
+        Repeating loading, start from negative slips
+
+        s : float
+        times : int
+        """
+        db = self.db
+        self.initialize_db(db)
+
+        # Monotonic pull-out
+        N = 0
+        smin = 0; smax = 0
+        slip = -self.slip_mono
+        bond = -self.bond_mono
+
+        self.store_to_db(db, slip, bond, smin, smax)
+
+        # Start the cycles
+        for i in range(0, times):
+            N += 1
+            t = self.get_bond(slip, bond, -s)
+            slip, bond, smin = self.define_bond_reloading(
+                slip, bond, -s, t, smin, smax, N)
+            self.store_to_db(db, slip, bond, smin, smax)
+
+            N += 1
+            t = self.get_bond(slip, bond, 0)
+            slip, bond, smax = self.define_bond_unloading(
+                slip, bond, 0, t, smin, smax, N)
+            self.store_to_db(db, slip, bond, smin, smax)
+
+        # End with a pull-out
+        N += 1
+        t = self.get_bond(slip, bond, -s)
+        slip, bond, smin = self.define_bond_reloading(
+            slip, bond, -s, t, smin, smax, N)
+        self.store_to_db(db, slip, bond, smin, smax)
+
+class TestEligenhausen(Eligenhausen):
+    """Tests Eligenhausen's bond-slip law"""
+    def __init__(self, fck, fctm):
+        """"""
+        super(TestEligenhausen, self).__init__(fck, fctm)
+
+        self.db = {}
+
+    def plot_bond(self, numbers, x_label="Slip (mm)", y_label = "Bond Stress (MPa)"):
+        """
+        Plots the local bond - local slip history for node i.
+        """
+        figure(0)
+        title("Eligenhausen's bond-slip law")
+        plot(self.slip_mono, self.bond_mono, label="Pos Mono")
+        plot(-self.slip_mono, -self.bond_mono, label="Neg Mono")
+
+        for j in numbers:
+            plot(self.db["slip_tot"][j], self.db["bond_tot"][j], label=str(j))
+        axhline(y=0, color = "black")
+        axvline(x=0, color = "black")
+        ylabel(y_label)
+        xlabel(x_label)
+        legend(loc = "best")
+        grid(True)
+
+        show()
+
+
+    def define_bond_monotonic(self):
+        self.smin = 0
+        self.tmin = 0
+        self.smax = 0
+        self.tmax = 0
+        self.E1 = 0
+        self.E2 = 0
+        self.Ef = 0
+        self.Esum = 0
+        self.Efsum = 0
+        self.d = 0
+        self.g = 0
+        self.df = 0
+        self.tf = 0
+        self.slip = self.slip_mono
+        self.bond = self.bond_mono
+
+    def recycling_tests_positive(self, s, times):
+        """
+        s : float
+        times : int
+        """
+        db = self.db
+        self.initialize_db(db)
+
+        # Monotonic pull-out
+        self.define_bond_monotonic()
+        self.store_to_db(self.db, self.slip, self.bond, self.smin, self.tmin,
+                         self.smax, self.tmax, self.E1, self.E2, self.Esum,
+                         self.d, self.g, self.Ef, self.Efsum, self.df, self.tf)
+
+        # Start the cycles
+        for i in range(0, times):
+            t = self.get_bond(self.slip, self.bond, s)
+
+            (self.smax, self.tmax, self.E1, self.E2, self.Esum, self.d, self.g,
+            self.df, self.tf, self.Ef, self.Efsum, self.slip, self.bond) = \
+                self.define_bond_unloading(s, t, self.slip, self.bond,
+                                           self.smin, self.smax, self.tmin,
+                                           self.Esum, self.Efsum)
+
+            self.store_to_db(db, self.slip, self.bond, self.smin, self.tmin,
+                             self.smax, self.tmax, self.E1, self.E2, self.Esum,
+                             self.d, self.g, self.Ef, self.Efsum, self.df, self.tf)
+
+            t = self.get_bond(self.slip, self.bond, -s)
+            (self.smin, self.tmin, self.E1, self.E2, self.Esum, self.d, self.g,
+            self.df, self.tf, self.Ef, self.Efsum, self.slip, self.bond) = \
+                self.define_bond_reloading(-s, t, self.slip, self.bond,
+                                           self.smin, self.smax, self.tmax,
+                                           self.Esum, self.Efsum)
+
+            self.store_to_db(db, self.slip, self.bond, self.smin, self.tmin,
+                             self.smax, self.tmax, self.E1, self.E2, self.Esum,
+                             self.d, self.g, self.Ef, self.Efsum, self.df, self.tf)
+        # finish with a pull-out
+        t = self.get_bond(self.slip, self.bond, s)
+
+        (self.smax, self.tmax, self.E1, self.E2, self.Esum, self.d, self.g,
+        self.df, self.tf, self.Ef, self.Efsum, self.slip, self.bond) = \
+            self.define_bond_unloading(s, t, self.slip, self.bond,
+                                       self.smin, self.smax, self.tmin,
+                                       self.Esum, self.Efsum)
+
+        self.store_to_db(db, self.slip, self.bond, self.smin, self.tmin,
+                         self.smax, self.tmax, self.E1, self.E2, self.Esum,
+                         self.d, self.g, self.Ef, self.Efsum, self.df, self.tf)
+    #def recycling_tests_negative(self, s, times):
+        #"""
+        #s : float
+        #times : int
+        #"""
+        #db = self.db
+        #self.initialize_db(db)
+
+        ## Monotonic pull-out
+        #smin = 0; tmin = 0; smax = 0; tmax = 0;
+        #E1 = 0; E2 = 0; Ef = 0; Esum = 0; Efsum = 0
+        #d = 0; g = 0; df = 0; tf = 0
+        #slip = -self.slip_mono
+        #bond = -self.bond_mono
+
+        #self.store_to_db(db, slip, bond, smin, tmin, smax, tmax,
+                           #E1, E2, Esum, d, g, df, tf, Ef, Efsum)
+
+        ## Start the cycles
+        #for i in range(0, times):
+            #t = self.get_bond(slip, bond, -s)
+            #slip, bond, smin, smax, tmin, Esum, Efsum = self.define_bond_reloading(
+                #-s, t, slip, bond, smin, smax, tmax, Esum, Efsum)
+
+            #t = self.get_bond(slip, bond, s)
+            #slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_unloading(
+                #s, t, slip, bond, smin, smax, tmin, Esum, Efsum)
+
+        ## end with a push-in
+        #t = self.get_bond(slip, bond, -s)
+        #slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_reloading(
+            #-s, t, slip, bond, smin, smax, tmin, Esum, Efsum)
+    #def repeating_tests_positive(self, s, times):
+        #"""
+        #s : float
+        #times : int
+        #"""
+        #db = self.db
+        #self.initialize_db(db)
+
+        ## Monotonic pull-out
+        #smin = 0; tmin = 0; smax = 0; tmax = 0;
+        #E1 = 0; E2 = 0; Ef = 0; Esum = 0; Efsum = 0
+        #d = 0; g = 0; df = 0; tf = 0
+        #slip = self.slip_mono
+        #bond = self.bond_mono
+
+        #self.store_to_db(db, slip, bond, smin, tmin, smax, tmax,
+                           #E1, E2, Esum, d, g, df, tf, Ef, Efsum)
+
+        ## Start the cycles
+        #for i in range(0, times):
+            #t = self.get_bond(slip, bond, s)
+            #slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_unloading(
+                #s, t, slip, bond, smin, smax, tmin, Esum, Efsum)
+
+            #t = self.get_bond(slip, bond, 0)
+            #slip, bond, smin, smax, tmin, Esum, Efsum = self.define_bond_reloading(
+                #0, t, slip, bond, smin, smax, tmax, Esum, Efsum)
+
+        ## end with a pull-out
+        #t = self.get_bond(slip, bond, s)
+        #slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_unloading(
+            #s, t, slip, bond, smin, smax, tmin, Esum, Efsum)
+
+    #def repeating_tests_negative(self, s, times):
+        #"""
+        #s : float
+        #times : int
+        #"""
+        #db = self.db
+        #self.initialize_db(db)
+
+        ## Monotonic pull-out
+        #smin = 0; tmin = 0; smax = 0; tmax = 0;
+        #E1 = 0; E2 = 0; Ef = 0; Esum = 0; Efsum = 0
+        #d = 0; g = 0; df = 0; tf = 0
+        #slip = -self.slip_mono
+        #bond = -self.bond_mono
+
+        #self.store_to_db(db, slip, bond, smin, tmin, smax, tmax,
+                           #E1, E2, Esum, d, g, df, tf, Ef, Efsum)
+
+        ## Start the cycles
+        #for i in range(0, times):
+            #t = self.get_bond(slip, bond, -s)
+            #slip, bond, smin, smax, tmin, Esum, Efsum = self.define_bond_reloading(
+                #-s, t, slip, bond, smin, smax, tmax, Esum, Efsum)
+
+            #t = self.get_bond(slip, bond, 0)
+            #slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_unloading(
+                #0, t, slip, bond, smin, smax, tmin, Esum, Efsum)
+
+        ## end with a push-in
+        #t = self.get_bond(slip, bond, -s)
+        #slip, bond, smin, smax, tmax, Esum, Efsum = self.define_bond_reloading(
+            #-s, t, slip, bond, smin, smax, tmin, Esum, Efsum)
+
+
+def test_bond_Tassios():
+    C30 = ModelCode2010(30)
+
+    s = 1
+
+    t = TestTassios(C30.fck, C30.fctm)
+    t.recycling_tests_positive(s, 10)
+    t.plot_bond((1, 2, 3, 4))
+    #t.recycling_tests_negative(s, 10)
+    #t.plot_bond((1, 2, 3, 4))
+    #t.repeating_tests_positive(s, 10)
+    #t.plot_bond((1, 2, 3, 4))
+    #t.repeating_tests_negative(s, 10)
+    #t.plot_bond((1, 2, 3, 4))
+    1
+
+def test_bond_Eligenhausen():
+    C30 = ModelCode2010(30)
+    T = TestEligenhausen(C30.fck, C30.fctm)
+
+    ss = (4.6,) #, 1.65, 4.6) #(0.009987, 0.01043, 0.0994, 0.1094)
+
+    for s in ss:
+        T.recycling_tests_positive(s, 10)
+        T.plot_bond((1,20))
+        #T.recycling_tests_negative(s, 10)
+        #T.plot_bond((1,2,3,4))
+        #T.repeating_tests_positive(s, 10)
+        #T.plot_bond((1,2,3,4))
+        #T.repeating_tests_negative(s, 10)
+        #T.plot_bond((1,2,3,4))
+    1
+
+if __name__ == "__main__":
+    test_bond_Tassios()
+    test_bond_Eligenhausen()
+    1

bunch.py

-#!/usr/bin/env python
-# encoding: utf-8
-"""
-Bunch is a subclass of dict with attribute-style access.
-
-from https://github.com/dsc/bunch/blob/master/bunch/__init__.py
-
->>> b = Bunch()
->>> b.hello = 'world'
->>> b.hello
-'world'
->>> b['hello'] += "!"
->>> b.hello
-'world!'
->>> b.foo = Bunch(lol=True)
->>> b.foo.lol
-True
->>> b.foo is b['foo']
-True
-
-It is safe to import * from this module:
-__all__ = ('Bunch', 'bunchify','unbunchify')
-un/bunchify provide dictionary conversion; Bunches can also be
-converted via Bunch.to/fromDict().
-"""
-
-__all__ = ('Bunch', 'bunchify','unbunchify')
-
-
-class Bunch(dict):
-    """
-    A dictionary that provides attribute-style access.
-    >>> b = Bunch()
-    >>> b.hello = 'world'
-    >>> b.hello
-    'world'
-    >>> b['hello'] += "!"
-    >>> b.hello
-    'world!'
-    >>> b.foo = Bunch(lol=True)
-    >>> b.foo.lol
-    True
-    >>> b.foo is b['foo']
-    True
-    
-    A Bunch is a subclass of dict; it supports all the methods a dict does...
-    >>> b.keys()
-    ['foo', 'hello']
-    
-    Including update()...
-    >>> b.update({ 'ponies': 'are pretty!' }, hello=42)
-    >>> print repr(b)
-    Bunch(foo=Bunch(lol=True), hello=42, ponies='are pretty!')
-    
-    As well as iteration...
-    >>> [ (k,b[k]) for k in b ]
-    [('ponies', 'are pretty!'), ('foo', Bunch(lol=True)), ('hello', 42)]
-    
-    And "splats".
-    >>> "The {knights} who say {ni}!".format(**Bunch(knights='lolcats', ni='can haz'))
-    'The lolcats who say can haz!'
-    
-    See unbunchify/Bunch.toDict, bunchify/Bunch.fromDict for notes about conversion.
-    """
-
-    def __contains__(self, k):
-        """
-        >>> b = Bunch(ponies='are pretty!')
-        >>> 'ponies' in b
-        True
-        >>> 'foo' in b
-        False
-        >>> b['foo'] = 42
-        >>> 'foo' in b
-        True
-        >>> b.hello = 'hai'
-        >>> 'hello' in b
-        True
-        """
-        try:
-            return hasattr(self, k) or dict.__contains__(self, k)
-        except:
-            return False
-
-    # only called if k not found in normal places
-    def __getattr__(self, k):
-        """
-        Gets key if it exists, otherwise throws AttributeError.
-        nb. __getattr__ is only called if key is not found in normal places.
-        >>> b = Bunch(bar='baz', lol={})
-        >>> b.foo
-        Traceback (most recent call last):
-        ...
-        AttributeError: foo
-        >>> b.bar
-        'baz'
-        >>> getattr(b, 'bar')
-        'baz'
-        >>> b['bar']
-        'baz'
-        >>> b.lol is b['lol']
-        True
-        >>> b.lol is getattr(b, 'lol')
-        True
-        """
-        try:
-            return self[k]
-        except KeyError:
-            raise AttributeError(k)
-
-    def __setattr__(self, k, v):
-        """
-        Sets attribute k if it exists, otherwise sets key k. A KeyError
-        raised by set-item (only likely if you subclass Bunch) will
-        propagate as an AttributeError instead.
-        >>> b = Bunch(foo='bar', this_is='useful when subclassing')
-        >>> b.values #doctest: +ELLIPSIS
-        <built-in method values of Bunch object at 0x...>
-        >>> b.values = 'uh oh'
-        >>> b.values
-        'uh oh'
-        >>> b['values']
-        Traceback (most recent call last):
-        ...
-        KeyError: 'values'
-        """
-        try:
-            # Throws exception if not in prototype chain
-            object.__getattribute__(self, k)
-        except AttributeError:
-            try:
-                self[k] = v
-            except:
-                raise AttributeError(k)
-        else:
-            object.__setattr__(self, k, v)
-
-    def __delattr__(self, k):
-        """
-        Deletes attribute k if it exists, otherwise deletes key k. A KeyError
-        raised by deleting the key--such as when the key is missing--will
-        propagate as an AttributeError instead.
-        >>> b = Bunch(lol=42)
-        >>> del b.values
-        Traceback (most recent call last):
-        ...
-        AttributeError: 'Bunch' object attribute 'values' is read-only
-        >>> del b.lol
-        >>> b.lol
-        Traceback (most recent call last):
-        ...
-        AttributeError: lol
-        """
-        try:
-            # Throws exception if not in prototype chain
-            object.__getattribute__(self, k)
-        except AttributeError:
-            try:
-                del self[k]
-            except KeyError:
-                raise AttributeError(k)
-        else:
-            object.__delattr__(self, k)
-
-    def toDict(self):
-        """
-        Recursively converts a bunch back into a dictionary.
-        >>> b = Bunch(foo=Bunch(lol=True), hello=42, ponies='are pretty!')
-        >>> b.toDict()
-        {'ponies': 'are pretty!', 'foo': {'lol': True}, 'hello': 42}
-        
-        See unbunchify for more info.
-        """
-        return unbunchify(self)
-
-    def __repr__(self):
-        """
-        Invertible* string-form of a Bunch.
-        >>> b = Bunch(foo=Bunch(lol=True), hello=42, ponies='are pretty!')
-        >>> print repr(b)
-        Bunch(foo=Bunch(lol=True), hello=42, ponies='are pretty!')
-        >>> eval(repr(b))
-        Bunch(foo=Bunch(lol=True), hello=42, ponies='are pretty!')
-        
-        (*) Invertible so long as collection contents are each repr-invertible.
-        """
-        keys = self.keys()
-        keys.sort()
-        args = ', '.join(['%s=%r' % (key, self[key]) for key in keys])
-        return '%s(%s)' % (self.__class__.__name__, args)
-
-    @staticmethod
-    def fromDict(d):
-        """
-        Recursively transforms a dictionary into a Bunch via copy.
-        >>> b = Bunch.fromDict({'urmom': {'sez': {'what': 'what'}}})
-        >>> b.urmom.sez.what
-        'what'
-        
-        See bunchify for more info.
-        """
-        return bunchify(d)
-
-
-
-# While we could convert abstract types like Mapping or Iterable, I think
-# bunchify is more likely to "do what you mean" if it is conservative about
-# casting (ex: isinstance(str,Iterable) == True ).
-#
-# Should you disagree, it is not difficult to duplicate this function with
-# more aggressive coercion to suit your own purposes.
-
-def bunchify(x):
-    """
-    Recursively transforms a dictionary into a Bunch via copy.
-    >>> b = bunchify({'urmom': {'sez': {'what': 'what'}}})
-    >>> b.urmom.sez.what
-    'what'
-    
-    bunchify can handle intermediary dicts, lists and tuples (as well as
-    their subclasses), but ymmv on custom datatypes.
-    >>> b = bunchify({ 'lol': ('cats', {'hah':'i win again'}), 'hello': [{'french':'salut', 'german':'hallo'}] })
-    >>> b.hello[0].french
-    'salut'
-    >>> b.lol[1].hah
-    'i win again'
-    
-    nb. As dicts are not hashable, they cannot be nested in sets/frozensets.
-    """
-    if isinstance(x, dict):
-        return Bunch( (k, bunchify(v)) for k,v in x.iteritems() )
-    elif isinstance(x, (list, tuple)):
-        return type(x)( bunchify(v) for v in x )
-    else:
-        return x
-
-def unbunchify(x):
-    """
-    Recursively converts a Bunch into a dictionary.
-    >>> b = Bunch(foo=Bunch(lol=True), hello=42, ponies='are pretty!')
-    >>> unbunchify(b)
-    {'ponies': 'are pretty!', 'foo': {'lol': True}, 'hello': 42}
-    
-    unbunchify will handle intermediary dicts, lists and tuples (as well as
-    their subclasses), but ymmv on custom datatypes.
-    >>> b = Bunch(foo=['bar', Bunch(lol=True)], hello=42, ponies=('are pretty!', Bunch(lies='are trouble!')))
-    >>> unbunchify(b)
-    {'ponies': ('are pretty!', {'lies': 'are trouble!'}), 'foo': ['bar', {'lol': True}], 'hello': 42}
-    
-    nb. As dicts are not hashable, they cannot be nested in sets/frozensets.
-    """
-    if isinstance(x, dict):
-        return dict( (k, unbunchify(v)) for k,v in x.iteritems() )
-    elif isinstance(x, (list, tuple)):
-        return type(x)( unbunchify(v) for v in x )
-    else:
-        return x
-
-
-if __name__ == "__main__":
-    import doctest
-    doctest.testmod()
 #! /usr/bin/env python2
 # -*- coding: utf-8 -*-
-# module : new.py
+# module : concrete.py
 ##------------------------------------------------------------------------------
 
 ##------------------------------------------------------------------------------
 # Licence : GPL v3
 ##------------------------------------------------------------------------------
 
+"""
+module : concrete.py
+
+Contains classes defining concrete's properties.
+"""
+
+##------------------------------------------------------------------------------
+# Importing necessary libraries
+##------------------------------------------------------------------------------
 from __future__ import division
 from __future__ import print_function
-from math import sqrt, log
-import numpy
-from matplotlib.pyplot import (figure, axhline, axvline, title, legend,
-                               show, xlabel, ylabel, grid, plot,
-                               close as closefig)
+from __future__ import unicode_literals
+
+from math import sqrt
+from math import log
 
 class Concrete(object):
     """An abstract concrete class"""
-    def __init__(self, fck, fctm = 3):
+    def __init__(self):
         super(Concrete, self).__init__()
 
-        self.fck = 20
-        self.fctm = fctm
-
-    def get_concrete_strain(self, Sc):
+    def get_strain(self, Sc):
         """Placeholder - must be overriden by a subclass."""
         raise NotImplementedError
 
         TODO
         Στο MC λαμβάνεται μειωμένο το E. σελίδα 119
         """
-        super(ModelCode2010, self).__init__(fck)
+        super(ModelCode2010, self).__init__()
         self.Df = Df
         self.vc = v
         self.ae = ae
 
         self.ecu = 0.002
 
-    def get_concrete_strain (self, Sc):
+
+    def get_strain (self, Sc):
         """
         Returns the concrete strain that corresponds to stress equal to Sc.
 
 
         return ec
 
-class TestConcrete(ModelCode2010):
-    def __init__(self, fck):
-        super(TestConcrete, self).__init__(fck)
-    def plot_concrete(self, x_label="ec (o/oo)", y_label="Sc (MPa)"):
-        Sc = numpy.linspace(self.fctm, -self.fck, 100)
-        ec = numpy.empty(len(Sc))
-        for i in range(len(Sc)):
-            ec[i] = self.get_concrete_strain(Sc[i])
-
-        figure(0)
-        title("Model Code 2010 strain - stress")
-        plot(-ec, -Sc)
-        axhline(y=0, color = "black")
-        axvline(x=0, color = "black")
-        ylabel(y_label)
-        xlabel(x_label)
-        grid(True)
-
-        show()
-
-def test():
-    C30 = TestConcrete(30)
-    C30.plot_concrete()
-
-if __name__ == "__main__":
-    test()
+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : concrete_test.py
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+"""
+module : concrete_test.py
+
+Contains test classes for concrete's models
+"""
+
+##------------------------------------------------------------------------------
+# Importing necessary libraries
+##------------------------------------------------------------------------------
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import unittest
+import numpy
+from matplotlib.pyplot import (figure, axhline, axvline, title, legend,
+                               show, xlabel, ylabel, grid, plot,
+                               close as closefig)
+
+from concrete import ModelCode2010
+
+class PlotConcrete(ModelCode2010):
+    def __init__(self, fck):
+        super(PlotConcrete, self).__init__(fck)
+
+    def plot_concrete(self, x_label="ec (o/oo)", y_label="Sc (MPa)"):
+        """
+        Plots concrete constitutive curve (σ - ε)
+        """
+        Sc = numpy.linspace(self.fctm, -self.fck, 100)
+        ec = numpy.empty(len(Sc))
+        for i in range(len(Sc)):
+            ec[i] = self.get_strain(Sc[i])
+
+        figure(0)
+        title("Model Code 2010 strain - stress")
+        plot(-ec, -Sc)
+        axhline(y=0, color = "black")
+        axvline(x=0, color = "black")
+        ylabel(y_label)
+        xlabel(x_label)
+        grid(True)
+
+        show()
+
+def test():
+    C30 = PlotConcrete(30.)
+    C30.plot_concrete()
+    print(C30.fck)
+    print(C30.Eci)
+
+if __name__ == "__main__":
+    test()

constants.py

-#! /usr/bin/env python
-# -*- coding: utf-8 -*-
-# module : constants.py
-##------------------------------------------------------------------------------
-
-##------------------------------------------------------------------------------
-# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
-# Licence : GPL v3
-##------------------------------------------------------------------------------
-
-"""
-constants.py Module
-
-This module contains the constants of the program.
-"""
-
-##------------------------------------------------------------------------------
-# Importing necessary libraries
-##------------------------------------------------------------------------------
-from __future__ import division
-
-##------------------------------------------------------------------------------
-# Defining constants
-##------------------------------------------------------------------------------
-## Number of points defining the bond-slip constitutive law
-BOND_POINTS = 10
-## The maximum number of iterations allowed in order to achieve convergence
-MAX_ITERATIONS = 500
-## The tolerance for which convergence is assumed.
-TOLERANCE = 1E-20
-## The slip value for which we suppose that there is significant slippage
-SIGNIFICANT_SLIP = 0.001
-## The length of the infinitestimal element in meters. 1/100 is one cm.
-## For development 1/100 is suggested. For final analysis 1/1000.
-INFINITESIMAL_LENGTH = 1/100
-## the length (in m) where we suppose an omoiothetic decrease of the bond-slip law.
-LENGTH_OF_OMOIOTHECY = 0.05
-## The name of the folder where the data of the analysis are going to be pickled
-DATA_FOLDER = "data"
-## If True then omoiothecy length is taken under consideration for the front of the specimen.
-OMOIOTHECY_FRONT = True
-#OMOIOTHECY_FRONT = False
-## If True then omoiothecy length is taken under consideration for the back of the specimen.
-#OMOIOTHECY_BACK = True
-OMOIOTHECY_BACK = False
-## There are numerical instabilities when unloading to zero, so we stop at +-
-## LOAD_STABILITY.
-LOAD_STABILITY = 0.1
-##------------------------------------------------------------------------------
+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : lines.pxd
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+"""
+Cython declaration file.
+"""
+
+cimport cython
+cimport numpy
+from cpython cimport bool
+
+@cython.locals(pred = double,
+               item = double)
+cpdef bool isascending(list seq,
+                  bool strictly = *)
+
+@cython.locals(pred = double,
+               item = double)
+cpdef bool isdescending(list seq,
+                  bool strictly = *)
+
+@cython.locals(area = double,
+               i = int)
+cpdef double calc_area(list a,
+                       list b)
+
+@cython.locals(y0 = double,
+               i = int,
+               np = int)
+cpdef double interpolation(numpy.ndarray[double, ndim=1] xp,
+                           numpy.ndarray[double, ndim=1] yp,
+                           double x0)
-#! /usr/bin/env python
+#! /usr/bin/env python2
 # -*- coding: utf-8 -*-
 # module : lines.py
 ##------------------------------------------------------------------------------
 
 ##------------------------------------------------------------------------------
-# Author : P. Mavrogiorgos (pmav99@gmail.com)
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
 # Licence : GPL v3
 ##------------------------------------------------------------------------------
 
 # Importing necessary libraries
 ##------------------------------------------------------------------------------
 from __future__ import division
-import operator
+from __future__ import print_function
+from __future__ import unicode_literals
 
-def calc_area(a, b):
-    """
-    Calculates the area of an irregular convex polygon. The points must be
-    placed clockwise or counterclockwise. No convexivity checks are being done.
-
-    Notes
-    -----
-    When the points are clockwise, then the calculated area is negative, so the
-    absolute value of area is returned.
-    For a N-polygon you input N points. The first point doesn't need to be
-    inputed twice.
-
-    References
-    ----------
-    http://paulbourke.net/geometry/polyarea/
-    
-    Parameters
-    ----------
-    x : array like
-        The abscissas of the points
-    y : array like
-        The ordinates of the points
-
-    Returns
-    -------
-    area : float
-        The area of the polygon.
-
-    Test
-    ----
-    >>> import numpy as np
-
-    # Line segment
-    >>> a = np.array([0.00, 1.00])
-    >>> b = np.array([0.00, 1.00])
-    >>> calc_area(a, b)
-    0.0
-
-    # Triangle
-    >>> a = np.array([0.00, 1.00, 2.00])
-    >>> b = np.array([0.00, 1.00, 0.00])
-    >>> calc_area(a, b)
-    1.0
-
-    # 4 sides
-    >>> a = np.array([0.00, 1.00, 2.00, 3.00])
-    >>> b = np.array([0.00, 1.00, 1.00, 0.00])
-    >>> calc_area(a, b)
-    2.0
-
-    # Polygon
-    >>> x = np.array([0.72, 2.66, 5, 3.63, 4, 1.9, 1.9])
-    >>> y = np.array([2.28, 4.71, 3.5, 2.52, 1.6, 1, 1])
-    >>> calc_area(x, y)
-    8.3593000000000011
-
-    """
-    area = 0
-    for i in range(-1, len(a) - 1):
-        area += a[i] * b[i+1] - a[i+1] * b[i]
-    return abs(area) / 2
+from operator import le, ge, gt, lt
 
 def line_slope(xA, yA, xB, yB):
     """
         Returns `True` if the sequence's elements are in ascending order,
         `False` otherwise.
 
-    Test
-    ----
-    >>> import numpy as np
-    >>> x = np.array([-3, -2, 0, 2, 4])
-    >>> y = np.array([-3, -2, 0, 2, 2])
-    >>> isascending(x)
-    True
-    >>> isascending(y)
-    False
-    >>> isascending(y, strictly = False)
-    True
 
     """
     if strictly:
-        op = operator.ge
+        op = ge
     else:
-        op = operator.gt
+        op = gt
 
     pred = seq[0]
     for item in seq[1:]:
         Returns `True` if the sequence's elements are in descending order,
         `False` otherwise.
 
-    Test
-    ----
-    >>> import numpy as np
-    >>> x = np.array([4, 3, 2, 1, 0])
-    >>> y = np.array([4, 3, 2, 1, 1])
-    >>> isdescending(x)
-    True
-    >>> isdescending(y)
-    False
-    >>> isdescending(y, strictly = False)
-    True
-
     """
     if strictly:
-        op = operator.le
+        op = le
     else:
-        op = operator.lt
+        op = lt
 
     pred = seq[0]
     for item in seq[1:]:
     function. E.g. interpolation(yp, xp, y0) where y and x are lists
     containing the coordinates of the points defining the curve.
 
+    Parameters
+    ----------
+    xp : numpy array (1D)
+        the abscissas of the points describing the curve.
+    yp : numpy array (1D)
+        the ordinates of the points describing the curve.
+    x0 : double
+        the given abscissa.
 
-    * x0     : the given abscissa.
-    * xp : the abscissas of the points describing the curve.
-    * yp : the ordinates of the points describing the curve.
-
-    Output
-
-    Tests
-    -----
-    >>> import numpy as np
-
-    >>> asc_x = np.array([0, 1, 3, 5, 10])
-    >>> des_x = np.array([10, 5, 3, 1, 0])
-    >>> asc_y = np.array([2, 4, 6, 8, 10])
-    >>> des_y = np.array([10, 8, 6, 4, 2])
-
-    >>> test_values = [-2, 1, 3, 4.65, 7.5, 11]
-
-    >>> [interpolation(asc_x, asc_y, value) for value in test_values]
-    [-2.0, 4.0, 6.0, 7.6500000000000004, 9.0, 10.4]
-    >>> [interpolation(des_x, asc_y, value) for value in test_values]
-    [14.0, 8.0, 6.0, 4.3499999999999996, 3.0, 1.6000000000000001]
-    >>> [interpolation(asc_x, des_y, value) for value in test_values]
-    [14.0, 8.0, 6.0, 4.3499999999999996, 3.0, 1.6000000000000001]
-    >>> [interpolation(des_x, des_y, value) for value in test_values]
-    [-2.0, 4.0, 6.0, 7.6500000000000004, 9.0, 10.4]
+    Returns
+    -------
+    area : float
+        The area of the polygon.
 
     """
     # Calculate the number of points that define the curve.
     np = len(xp)
-    y0 = 999999999.
+    y0 = 999999999
 
     # check if the points are in ascending order.
     if xp[0] < xp[1]:
         ## If x0 > x_max interpolate between the last two points
         if x0 > xp[np-1]:
-            y0 = (x0 - xp[np-2]) *\
-               (yp[np-1] - yp[np-2]) /\
-               (xp[np-1] - xp[np-2]) + yp[np-2]
+            y0 = ((x0 - xp[np-2]) * (yp[np-1] - yp[np-2]) /
+                  (xp[np-1] - xp[np-2]) +
+                  yp[np-2])
         else:
             for i in range(1, np, 1):
                 if x0 <= xp[i]:
-                    y0 = (x0 - xp[i-1]) *\
-                       (yp[i] - yp[i-1]) /\
-                       (xp[i] - xp[i-1]) + yp[i-1]
+                    y0 = ((x0 - xp[i-1]) * (yp[i] - yp[i-1]) /
+                          (xp[i] - xp[i-1]) +
+                          yp[i-1])
                     break
 
     # check if the points are in descending order.
     elif xp[0] > xp[1]:
         # If x0 < x_min interpolate between the last two points
         if x0 < xp[np-1]:
-            y0 = (x0 - xp[np-2]) *\
-               (yp[np-1] - yp[np-2]) /\
-               (xp[np-1] - xp[np-2]) + yp[np-2]
+            y0 = ((x0 - xp[np-2]) * (yp[np-1] - yp[np-2]) /
+                  (xp[np-1] - xp[np-2]) +
+                  yp[np-2])
         else:
             for i in range(1 ,np, 1):
                 if x0 >= xp[i]:
-                    y0 = (x0 - xp[i-1]) *\
-                       (yp[i] - yp[i-1]) /\
-                       (xp[i] - xp[i-1]) + yp[i-1]
+                    y0 = ((x0 - xp[i-1]) * (yp[i] - yp[i-1]) /
+                          (xp[i] - xp[i-1]) +
+                          yp[i-1])
                     break
     else:
-        print "Wrong Point Values!!!"
+        print("Wrong Point Values!!!")
 
     return y0
 
-def calc_E1(x, y, x0, y0):
+def calc_area(a, b):
     """
-    >>> x = [-1, 0, 1, 2, 4]
-    >>> y = [-1, 0, 1, 1, 3]
-    >>> calc_E1(x, y, -5, 3)
-    0.0
-    >>> calc_E1(x, y, 0, 3)
-    0.0
-    >>> calc_E1(x, y, 1, 3)
-    0.5
-    >>> calc_E1(x, y, 3, 3)
-    3.5
-    >>> calc_E1(x, y, 6, 3)
-    Traceback (most recent call last):
-    ...
-    IOError: x0 out of bounds
+    Calculates the area of an irregular convex polygon. The points must be
+    placed clockwise or counterclockwise. No convexivity checks are being done.
+
+    Notes
+    -----
+    When the points are clockwise, then the calculated area is negative, so the
+    absolute value of area is returned.
+    For a N-polygon you input N points. The first point doesn't need to be
+    inputed twice.
+
+    References
+    ----------
+    http://paulbourke.net/geometry/polyarea/
+
+    Parameters
+    ----------
+    x : list
+        The abscissas of the points
+    y : list
+        The ordinates of the points
+
+    Returns
+    -------
+    area : float
+        The area of the polygon.
+
     """
-    # slip value must be lower than the max slip curve value
-    if x0 > x[-1]:
-        raise IOError("x0 out of bounds")
-
-    # bond's curve second point must be 0!
-    if y[1] != 0:
-        raise IOError("wrong values in curve definition.")
-
     area = 0
-    if y[0] < 0:
-        for j, value in enumerate(x):
-            if x0 < value:
-                break
-
-        a = list(x[1:j])
-        b = list(y[1:j])
-        a.append(x0)
-        b.append(y0)
-        a.append(x0)
-        b.append(0)
-
-        for i in range(-1, len(a) - 1):
-            area += a[i] * b[i+1] - a[i+1] * b[i]
-
-    return abs(area) / 2
-
-if __name__ == "__main__":
-    import doctest
-    doctest.testmod(verbose = False)
+    for i in range(-1, len(a) - 1):
+        area += a[i] * b[i+1] - a[i+1] * b[i]
+    return abs(area) / 2
+import sys
+
+# setup.py
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Distutils import build_ext
+import numpy
+
+# make setup build extension modules in place
+sys.argv.append("build_ext")
+sys.argv.append("--inplace")
+
+# for notes on compiler flags see:
+# http://docs.python.org/install/index.html
+
+cmdclass = {"build_ext": build_ext}
+
+include_dirs = [numpy.get_include()]
+
+ext_modules = [Extension("lines",
+                         ["lines.py"],
+                         extra_compile_args=["-O3"],
+                         )
+               ]
+
+for e in ext_modules:
+    e.pyrex_directives = {"embedsignature": True}
+
+setup(
+    cmdclass = cmdclass,
+    include_dirs = include_dirs,
+    ext_modules = ext_modules
+)
+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : lines_test.py
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+"""
+lines_test.py module
+
+Contains unittests for the functions defined in lines.py
+(isascending, isdescending, interpolation, calc_area)
+"""
+
+##------------------------------------------------------------------------------
+# Importing necessary libraries
+##------------------------------------------------------------------------------
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import unittest
+from numpy import array
+
+from lines import isascending
+from lines import isdescending
+from lines import interpolation
+from lines import calc_area
+
+class IsAscendingTest(unittest.TestCase):
+    def setUp(self):
+        self.asc_strict = [-3, -2, 0, 2, 4]
+        self.asc = [-3, -2, 0, 2, 2]
+        self.desc = [10, 9, 8, 7, 5]
+
+    def test_isascendig_strictly_true(self):
+        self.assertTrue(isascending(self.asc_strict, strictly = True))
+
+    def test_isascendig_strictly_false(self):
+        self.assertFalse(isascending(self.desc, strictly = True))
+
+    def test_isascending_true(self):
+        self.assertTrue(isascending(self.asc, strictly = False))
+
+    def test_isascending_false(self):
+        self.assertFalse(isascending(self.desc, strictly = False))
+
+class IsDescendingTest(unittest.TestCase):
+    def setUp(self):
+        self.desc_strict = [4, 3, 2, 1, 0]
+        self.desc = [4, 3, 2, 1, 1]
+        self.asc = [0, 1, 2, 3, 4]
+
+    def test_isdescendig_strictly_true(self):
+        self.assertTrue(isdescending(self.desc_strict, strictly = True))
+
+    def test_isdescendig_strictly_false(self):
+        self.assertFalse(isdescending(self.desc, strictly = True))
+
+    def test_isdescending_true(self):
+        self.assertTrue(isdescending(self.desc, strictly = False))
+
+    def test_isdescending_false(self):
+        self.assertFalse(isdescending(self.asc, strictly = False))
+
+class InterpolationTest(unittest.TestCase):
+    def setUp(self):
+        self.asc_x = array([0., 1., 3., 5., 10.])
+        self.des_x = array([10., 5., 3., 1., 0.])
+        self.asc_y = array([2., 4., 6., 8., 10.])
+        self.des_y = array([10., 8., 6., 4., 2.])
+        #self.test_values = [-2, 1, 3, 4.65, 7.5, 11]
+
+    def test_asc_vs_des_lower(self):
+        self.assertEquals(14., interpolation(self.asc_x, self.des_y, -2.))
+
+    def test_asc_vs_des_within1(self):
+        self.assertEquals(8., interpolation(self.asc_x, self.des_y, 1.))
+
+    def test_asc_vs_des_within2(self):
+        self.assertEquals(6., interpolation(self.asc_x, self.des_y, 3.))
+
+    def test_asc_vs_des_within3(self):
+        self.assertAlmostEquals(4.34999999, interpolation(self.asc_x, self.des_y, 4.65))
+
+    def test_asc_vs_des_within4(self):
+        self.assertEquals(3., interpolation(self.asc_x, self.des_y, 7.5))
+
+    def test_asc_vs_des_upper(self):
+        self.assertEquals(1.6, interpolation(self.asc_x, self.des_y, 11.))
+
+    #############
+    def test_des_vs_des_lower(self):
+        self.assertEquals(-2., interpolation(self.des_x, self.des_y, -2))
+
+    def test_des_vs_des_within1(self):
+        self.assertEquals(4., interpolation(self.des_x, self.des_y, 1.))
+
+    def test_des_vs_des_within2(self):
+        self.assertEquals(6., interpolation(self.des_x, self.des_y, 3.))
+
+    def test_des_vs_des_within3(self):
+        self.assertAlmostEquals(7.65, interpolation(self.des_x, self.des_y, 4.65))
+
+    def test_des_vs_des_within4(self):
+        self.assertEquals(9., interpolation(self.des_x, self.des_y, 7.5))
+
+    def test_des_vs_des_upper(self):
+        self.assertEquals(10.4, interpolation(self.des_x, self.des_y, 11.))
+
+    #############
+    def test_des_vs_asc_lower(self):
+        self.assertEquals(14., interpolation(self.des_x, self.asc_y, -2.))
+
+    def test_des_vs_asc_within1(self):
+        self.assertEquals(8., interpolation(self.des_x, self.asc_y, 1.))
+
+    def test_des_vs_asc_within2(self):
+        self.assertEquals(6, interpolation(self.des_x, self.asc_y, 3.))
+
+    def test_des_vs_asc_within3(self):
+        self.assertAlmostEquals(4.34999999, interpolation(self.des_x, self.asc_y, 4.65))
+
+    def test_des_vs_asc_within4(self):
+        self.assertEquals(3., interpolation(self.des_x, self.asc_y, 7.5))
+
+    def test_des_vs_asc_upper(self):
+        self.assertEquals(1.6, interpolation(self.des_x, self.asc_y, 11.))
+
+    #############
+    def test_asc_vs_asc_lower(self):
+        self.assertEquals(-2., interpolation(self.asc_x, self.asc_y, -2.))
+
+    def test_asc_vs_asc_within1(self):
+        self.assertEquals(4., interpolation(self.asc_x, self.asc_y, 1.))
+
+    def test_asc_vs_asc_within2(self):
+        self.assertEquals(6., interpolation(self.asc_x, self.asc_y, 3.))
+
+    def test_asc_vs_asc_within3(self):
+        self.assertEquals(7.65, interpolation(self.asc_x, self.asc_y, 4.65))
+
+    def test_asc_vs_asc_within4(self):
+        self.assertEquals(9., interpolation(self.asc_x, self.asc_y, 7.5))
+
+    def test_asc_vs_asc_upper(self):
+        self.assertEquals(10.4, interpolation(self.asc_x, self.asc_y, 11.))
+
+class CalcAreaTest(unittest.TestCase):
+    def test_polygon1(self):
+        x = [0.72, 2.66, 5, 3.63, 4, 1.9, 1.9]
+        y = [2.28, 4.71, 3.5, 2.52, 1.6, 1, 1]
+        self.assertAlmostEquals(8.3593, calc_area(x, y))
+
+    def test_polygon2(self):
+        x = [1., 2, 3, 4, 4]
+        y = [2., 3, 4, 4, 2]
+        self.assertEquals(4., calc_area(x, y))
+
+    def test_4_sides(self):
+        x = [0.00, 1.00, 2.00, 3.00]
+        y = [0.00, 1.00, 1.00, 0.00]
+        self.assertEquals(2, calc_area(x, y))
+
+    def test_line_segment(self):
+        x = [0.00, 1.00]
+        y = [0.00, 1.00]
+        self.assertEquals(0, calc_area(x, y))
+
+    def test_triangle(self):
+        x = [0.00, 1.00, 2.00]
+        y = [0.00, 1.00, 0.00]
+        self.assertEquals(1, calc_area(x, y))
+
+def test_unittest():