Commits

Panagiotis Mavrogiorgos  committed 2cf247f

Major

  • Participants
  • Parent commits 17d9b37

Comments (0)

Files changed (14)

File Constants.py

 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/100
+INFINITESIMAL_LENGTH = 1/1000
 ## 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.

File bond_test.py

File contents unchanged.

File concrete.pxd

+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : concrete.pxd
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+"""
+Cython declaration file.
+"""
+
+cimport cython
+#cimport numpy
+#from cpython cimport bool
+
+#cdef extern from "math.h":
+    #double sqrt(double theta)
+    #double log(double theta)
+
+cdef class Concrete(object):
+    pass
+
+cdef class ModelCode2010(Concrete):
+    cpdef double fck, fctm, ecu, Eci
+
+    cpdef double get_strain (self, double Sc)

File concrete_setup.py

+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("concrete",
+                         ["concrete.py"],
+                         extra_compile_args=["-O3"],
+                         libraries=["m"],
+                         )
+               ]
+
+for e in ext_modules:
+    e.pyrex_directives = {"embedsignature": True}
+
+setup(
+    cmdclass = cmdclass,
+    include_dirs = include_dirs,
+    ext_modules = ext_modules
+)

File concrete_test.py

File contents unchanged.
+#!/usr/bin/env python
+#coding:utf-8
+# Author: Panos Mavrogiorgos -- pmav99 <gmail>
+# Purpose: Contours
+# Created: 10/18/2011
+
+from matplotlib.pyplot import (imshow, xlabel, ylabel, show, colorbar,clabel,
+                               figure, contour, contourf,grid)
+from matplotlib import rcParams
+from numpy import linspace, meshgrid
+
+# Set Up the matplotlib figures.
+rcParams['lines.linewidth'] = 1
+rcParams['font.size'] = 14
+rcParams['font.family'] = 'serif'
+rcParams['mathtext.fontset'] = 'stix'
+rcParams['axes.labelsize'] = 'x-large'
+rcParams['xtick.labelsize'] = 'medium'
+rcParams['ytick.labelsize'] = 'medium'
+rcParams['legend.fontsize'] = 'large'
+
+def plane_Z(X, Y):
+    temp = 0.0
+    a = 8.2263888888889547E-01
+    b = 5.7166666666666721E+00
+    c = 1.2875000000000014E-01
+    temp += a
+    temp += b * X
+    temp += c * Y
+    return temp
+
+def surface_Z(X, Y):
+    temp = 0.0
+    a = 1.2051388888888865E+00
+    b = -1.9333333333333216E+00
+    c = 8.6250000000000104E-02
+    d = 8.4999999999999920E-01
+
+    temp += a
+    temp += b * X
+    temp += c * Y
+    temp += d * X * Y
+    return temp
+
+def plane_Z_eligenhausen(X, Y):
+    temp = 0.0
+    a = 9.7527777777777813E-01
+    b = 1.3000000000000016E+00
+    c = 4.1666666666666519E-03
+
+    temp += a
+    temp += b * X
+    temp += c * Y
+    return temp
+
+def surface_Z_eligenhausen(X, Y):
+    temp = 0.0
+    a = 1.0090277777777767E+00
+    b = 6.2500000000000666E-01
+    c = 4.1666666666663639E-04
+    d = 7.4999999999999428E-02
+    temp += a
+    temp += b * X
+    temp += c * Y
+    temp += d * X * Y
+    return temp
+
+def plane(N = 100):
+    x = linspace(0, 0.10, N)
+    y = linspace(5, 13, N)
+    X, Y = meshgrid(x, y)
+    Z = plane_Z_eligenhausen(X, Y)
+
+    f = figure()
+    c = contour(X, Y, Z, colors = 'k')
+    clabel(c, fontsize = 12, inline = True)
+    xlabel("$q$")
+    ylabel("$N$")
+    grid(True)
+
+def surface(N = 100):
+    x = linspace(0, 0.10, N)
+    y = linspace(5, 13, N)
+    X, Y = meshgrid(x, y)
+    Z = surface_Z_eligenhausen(X, Y)
+
+    f = figure()
+    c = contour(X, Y, Z, colors = 'k')
+    clabel(c, fontsize = 12, inline = True)
+    xlabel("$q$")
+    ylabel("$N$")
+    grid(True)
+
+
+if __name__ == "__main__":
+    plane()
+    surface()
+    show()

File diereunisi.py

+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : new.py
+
+from __future__ import division
+from pull_out_tests import *
+from plot import *
+
+import shelve
+import cPickle
+
+"""
+Defaults
+--------
+INFINITESIMAL_LENGTH = 1/100
+LENGTH_OF_OMOIOTHECY = 0.05
+OMOIOTHECY_FRONT = True
+OMOIOTHECY_BACK = False
+"""
+
+#PullOutExternalTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc, lext = lext)
+#PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+
+import os
+import os.path
+import codecs
+
+from constants import *
+
+def tassios_Ds():
+    fck = 30
+    fyk = 500
+    P = 500
+    L = 1000
+    Dc = 120
+    #Ds = 20
+    #lext = 400
+
+    cwd = os.path.abspath(os.curdir)
+    fname = "Steel_Diameter"
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+    Ds12 = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = 12, Dc = Dc)
+    Ds12.run_recycling_positive(P, 1)
+    Ds14 = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = 14, Dc = Dc)
+    Ds14.run_recycling_positive(P, 1)
+    Ds16 = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = 16, Dc = Dc)
+    Ds16.run_recycling_positive(P, 1)
+    Ds18 = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = 18, Dc = Dc)
+    Ds18.run_recycling_positive(P, 1)
+    Ds20 = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = 20, Dc = Dc)
+    Ds20.run_recycling_positive(P, 1)
+    Ds25 = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = 25, Dc = Dc)
+    Ds25.run_recycling_positive(P, 1)
+    Ds30 = PullOutClassicTassios(fck = fck, fyk = fyk, L = 1400, Ds = 35, Dc = Dc)
+    Ds30.run_recycling_positive(P, 1)
+
+    latex = u""
+    s = u""
+    s += "fck = %d MPa\n" % fck
+    s += "fyk = %d MPa\n" % fyk
+    s += "P = %d MPa\n" % P
+    s += "L = %d mm\n" % L
+    s += "Dc = %d mm\n" % Dc
+    #s += "Ds = %d mm\n" % Ds
+    s += "Dx = %d mm\n" % int(INFINITESIMAL_LENGTH * 1000)
+    s += "Lom = %d mm \n" % int(LENGTH_OF_OMOIOTHECY * 1000)
+    s+= "-----------\n"
+    s+= "ds       n=1         n=3         s<0.001      s<0.0001\n"
+
+    # Το μήκος αγκύρωσης το εξετάζω μόνο για τη διαρροή.
+
+    for inst in (Ds12, Ds14, Ds16, Ds18, Ds20, Ds25, Ds30):
+        lb1 = 0
+        lb2 = 0
+        for i, slip in enumerate(inst.db.s_tot[1]):
+            if slip < 0.001 and lb1 == 0:
+                lb1 = int(i * INFINITESIMAL_LENGTH * 1000)
+            if slip < 0.0001:
+                lb2 = int(i * INFINITESIMAL_LENGTH * 1000)
+                break
+        s1 = inst.db.s_tot[1][0]
+        s3 = inst.db.s_tot[3][0]
+        s+= "%d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (int(inst.Ds*1000),
+                                                                               s1, s3, lb1, lb2,
+                                                                               s3 / s1)
+        latex += u"%d & %.3f & %.3f & %3d & %3d & %.2f\\\\\n" % (int(inst.Ds*1000),
+                                                                s1, s3, lb1, lb2,
+                                                                s3 / s1)
+
+    print s
+    print latex
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+    with codecs.open(os.path.join(fpath, fname + "_latex.txt"), "w", "utf-8") as f:
+        f.write(latex)
+
+    #if PLOT:
+        #plot.plot_slip_distribution("Ds monotonic",
+                               #(Ds12.db.s_tot[1], "Ds = 12 mm"),
+                               ##(Ds14.db.s_tot[1], "Ds = 14 mm"),
+                               #(Ds16.db.s_tot[1], "Ds = 16 mm"),
+                               ##(Ds18.db.s_tot[1], "Ds = 18 mm"),
+                               #(Ds20.db.s_tot[1], "Ds = 20 mm"))
+
+        #plot_slip_distribution("Ds After 1 cycle",
+                               #(Ds12.db.s_tot[3], "Ds = 12 mm, n = 3"),
+                               #(Ds14.db.s_tot[3], "Ds = 14 mm"),
+                               #(Ds16.db.s_tot[3], "Ds = 16 mm"),
+                               #(Ds18.db.s_tot[3], "Ds = 18 mm"),
+                               #(Ds20.db.s_tot[3], "Ds = 20 mm"))
+
+def tassios_steel():
+    fck = 30
+    #fyk = 500
+    #P = 500
+    L = 1000
+    Dc = 120
+    Ds = 20
+    #lext = 400
+
+    cwd = os.path.abspath(os.curdir)
+    fname = "Steel_strength"
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+    S400 = PullOutClassicTassios(fck = fck, fyk = 400, L = L, Ds = Ds, Dc = Dc)
+    S400.run_recycling_positive(400, 1)
+    S500 = PullOutClassicTassios(fck = fck, fyk = 500, L = L, Ds = Ds, Dc = Dc)
+    S500.run_recycling_positive(500, 1)
+    S600 = PullOutClassicTassios(fck = fck, fyk = 600, L = L, Ds = Ds, Dc = Dc)
+    S600.run_recycling_positive(600, 1)
+
+    latex = u""
+    s = u""
+    s += "fck = %d MPa\n" % fck
+    #s += "fyk = %d MPa\n" % fyk
+    #s += "P = %d MPa\n" % P
+    s += "L = %d mm\n" % L
+    s += "Dc = %d mm\n" % Dc
+    s += "Ds = %d mm\n" % Ds
+    s += "Dx = %d mm\n" % int(INFINITESIMAL_LENGTH * 1000)
+    s += "Lom = %d mm \n" % int(LENGTH_OF_OMOIOTHECY * 1000)
+    s+= "-----------\n"
+    s+= "fyk       n=1         n=3         s<0.001      s<0.0001\n"
+
+    # Το μήκος αγκύρωσης το εξετάζω μόνο για τη διαρροή.
+    for inst in (S400, S500, S600):
+        lb1 = 0
+        lb2 = 0
+        for i, slip in enumerate(inst.db.s_tot[1]):
+            if slip < 0.001 and lb1 == 0:
+                lb1 = int(i * INFINITESIMAL_LENGTH * 1000)
+            if slip < 0.0001:
+                lb2 = int(i * INFINITESIMAL_LENGTH * 1000)
+                break
+        s1 = inst.db.s_tot[1][0]
+        s3 = inst.db.s_tot[3][0]
+        s+= "%d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (inst.S.fy,
+                                                                               s1, s3, lb1, lb2,
+                                                                               s3 / s1)
+        latex += u"%d & %.3f & %.3f & %3d & %3d & %.2f\\\\\n" % (inst.S.fy,
+                                                                s1, s3, lb1, lb2,
+                                                                s3 / s1)
+
+    print s
+    print latex
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+    with codecs.open(os.path.join(fpath, fname + "_latex.txt"), "w", "utf-8") as f:
+        f.write(latex)
+
+def tassios_infinitesimal_length():
+    P = 500
+    fck = 30
+    fyk = 500
+    Ds = 20
+    Dc = Ds * 6
+    lext = 400
+    L = 1000
+
+    cwd = os.path.abspath(os.curdir)
+    fname = "Infinitesimal_length"
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+    latex = u""
+    s = u""
+    s+= "fyk = %d\nfck = %d\nDs = %d\nDc = %d\nP = %d\nL = %d\nLom = %d\n" % (fyk, fck, Ds, Dc, P, L, int(LENGTH_OF_OMOIOTHECY * 1000))
+    s+= "-----------\n"
+    s+= "inf       n=1         n=3         s<0.001      s<0.0001\n"
+
+    for inf in (1/10, 1/20, 1/50, 1/100, 1/200, 1/500, 1/1000, 1/2000):
+
+        inst = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+        inst.run_recycling_positive(500, 1)
+
+        lb1 = 0
+        lb2 = 0
+        for i, slip in enumerate(inst.db.s_tot[1]):
+            if slip < 0.001 and lb1 == 0:
+                lb1 = int(i * inst.Dx * 1000)
+            if slip < 0.0001:
+                lb2 = int(i * inst.Dx * 1000)
+                break
+        s1 = inst.db.s_tot[1][0]
+        s3 = inst.db.s_tot[3][0]
+        s+= "%d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (int(inst.Dx * 1000),
+                                                                                s1, s3, lb1, lb2,
+                                                                                s3 / s1)
+        latex += u"%.1f & %.3f & %.3f & %3d & %3d & %.2f\\\\\n" % (inst.Dx * 1000,
+                                                                 s1, s3, lb1, lb2,
+                                                             s3 / s1)
+    print s
+    print latex
+
+    #plot_omoiothecy_length(fpath)
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+    with codecs.open(os.path.join(fpath, fname + "_latex.txt"), "w", "utf-8") as f:
+        f.write(latex)
+
+    #def lom():
+        #for inf in (1/10, 1/20, 1/50, 1/100, 1/200, 1/500, 1/1000, 1/2000):
+            #yield inf
+
+
+    #a = lom()
+
+    # These must be added at Data __init__ method
+        #global LENGTH_OF_OMOIOTHECY, OMOIOTHECY_FRONT
+        #LENGTH_OF_OMOIOTHECY = a.next()
+        #if LENGTH_OF_OMOIOTHECY == 0.00001:
+            #OMOIOTHECY_FRONT = False
+        #else:
+            #OMOIOTHECY_FRONT = True
+
+        #global INFINITESIMAL_LENGTH
+        #INFINITESIMAL_LENGTH = a.next()
+
+def tassios_length():
+    L800_S500_C30_ds20_Dc120 = PullOutClassicTassios(
+        fck = 30, fyk = 500, L = 800, Ds = 20, Dc = 120)
+    L800_S500_C30_ds20_Dc120.run_recycling_positive(500, 1)
+    t1 = L800_S500_C30_ds20_Dc120
+
+    L1000_S500_C30_ds20_Dc120 = PullOutClassicTassios(
+        fck = 30, fyk = 500, L = 1000, Ds = 20, Dc = 120)
+    L1000_S500_C30_ds20_Dc120.run_recycling_positive(500, 1)
+    t2 = L1000_S500_C30_ds20_Dc120
+
+    L1200_S500_C30_ds20_Dc120 = PullOutClassicTassios(
+        fck = 30, fyk = 500, L = 1200, Ds = 20, Dc = 120)
+    L1200_S500_C30_ds20_Dc120.run_recycling_positive(500, 1)
+    t3 = L1200_S500_C30_ds20_Dc120
+
+    print(t1.db.s_tot[1][0], t1.db.s_tot[3][0])
+    print(t2.db.s_tot[1][0], t2.db.s_tot[3][0])
+    print(t3.db.s_tot[1][0], t3.db.s_tot[3][0])
+
+    """
+    INFINITESIMAL_LENGTH = 1/100
+    ---------------------------------
+    L700  : Not solved
+    L800  : 0.603885683656 0.779008736379 -------- k2 = 1.29
+    L1000 : 0.603885679852 0.779140238709 -------- k2 = 1.29
+    L1200 : 0.603885679849 0.779145889485 -------- k2 = 1.29
+
+    INFINITESIMAL_LENGTH = 1/1000
+    ---------------------------------
+    L700  : Not solved
+    L800  : 0.594818604402 0.769283181106 -------- k2 = 1.29
+    L1000 : 0.594818600939 0.769400513063 -------- k2 = 1.29
+    L1200 : 0.594818600937 0.769405129653 -------- k2 = 1.29
+    """
+
+def tassios_omoiothecy_length_front():
+    P = 500
+    fck = 30
+    fyk = 500
+    Ds = 20
+    Dc = Ds * 6
+    lext = 400
+    L = 1000
+
+    cwd = os.path.abspath(os.curdir)
+    fname = "Omoiothecy_Front"
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+    inst = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    inst.run_recycling_positive(500, 1)
+    s = u""
+    s+= "fyk = %d\nfck = %d\nDs = %d\nDc = %d\nP = %d\nL = %d\n" % (fyk, fck, Ds, Dc, P, int(inst.L * 1000))
+    s+= "-----------\n"
+    s+= "Lom       n=1         n=3         s<0.001      s<0.0001\n"
+
+    # Το μήκος αγκύρωσης το εξετάζω μόνο για τη διαρροή.
+
+    lb1 = 0
+    lb2 = 0
+    for i, slip in enumerate(inst.db.s_tot[1]):
+        if slip < 0.001 and lb1 == 0:
+            lb1 = int(i * INFINITESIMAL_LENGTH * 1000)
+        if slip < 0.0001:
+            lb2 = int(i * INFINITESIMAL_LENGTH * 1000)
+            break
+    s1 = inst.db.s_tot[1][0]
+    s3 = inst.db.s_tot[3][0]
+    s+= "%2d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (0,
+                                                                            s1, s3, lb1, lb2,
+                                                                            s3 / s1)
+    latex = u"%4d & %.3f & %.3f & %3d & %3d & %.2f\\\\\n" % (0,
+                                                             s1, s3, lb1, lb2,
+                                                             s3 / s1)
+
+    # The pull_out_test.so must be removed first!!!
+    # The *.py must be used after adding the following code.
+    # These must be added at pull_out test.
+    # First we create an iterator
+
+    #def lom():
+        #for length in (0.00001, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10):
+            #yield length
+
+    #a = lom()
+
+    # These must be added at Data __init__ method
+        #global LENGTH_OF_OMOIOTHECY, OMOIOTHECY_FRONT
+        #LENGTH_OF_OMOIOTHECY = a.next()
+        #if LENGTH_OF_OMOIOTHECY == 0.00001:
+            #OMOIOTHECY_FRONT = False
+        #else:
+            #OMOIOTHECY_FRONT = True
+
+    for length in (0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10):
+        inst = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+        inst.run_recycling_positive(500, 1)
+
+        # Το μήκος αγκύρωσης το εξετάζω μόνο για τη διαρροή.
+
+        lb1 = 0
+        lb2 = 0
+        for i, slip in enumerate(inst.db.s_tot[1]):
+            if slip < 0.001 and lb1 == 0:
+                lb1 = int(i * INFINITESIMAL_LENGTH * 1000)
+            if slip < 0.0001:
+                lb2 = int(i * INFINITESIMAL_LENGTH * 1000)
+                break
+        s1 = inst.db.s_tot[1][0]
+        s3 = inst.db.s_tot[3][0]
+        s+= "%2d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (inst.L_om * INFINITESIMAL_LENGTH * 1000,
+                                                                                s1, s3, lb1, lb2,
+                                                                                s3 / s1)
+        latex += u"%4d & %.3f & %.3f & %3d & %3d & %.2f\\\\\n" % (inst.L_om * INFINITESIMAL_LENGTH * 1000,
+                                                                  s1, s3, lb1, lb2,
+                                                                  s3 / s1)
+    print s
+    print latex
+
+    plot_omoiothecy_length(fpath)
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+    with codecs.open(os.path.join(fpath, fname + "_latex.txt"), "w", "utf-8") as f:
+        f.write(latex)
+
+def compare_with_fortran():
+    L = 1000
+    fck = 30
+    fyk = 420
+    Ds = 12
+    Dc = 113
+    lext = 400
+
+    cwd = os.path.abspath(os.curdir)
+    #fname = "Classic_Ds%d" % Ds
+    fname = "External_Ds%d_Lb%d" % (Ds, lext)
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+    t420 = PullOutExternalTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc, lext = lext)
+    t420.run_recycling_positive(420, 1)
+
+    t300 = PullOutExternalTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc, lext = lext)
+    t300.run_recycling_positive(300, 1)
+
+    t200 = PullOutExternalTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc, lext = lext)
+    t200.run_recycling_positive(200, 1)
+
+    t100 = PullOutExternalTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc, lext = lext)
+    t100.run_recycling_positive(100, 1)
+
+    # Το μήκος αγκύρωσης το εξετάζω μόνο για τη διαρροή.
+    lb1 = 0
+    lb2 = 0
+    for i, s in enumerate(t420.db.s_tot[1]):
+        if s < 0.001 and lb1 == 0:
+            lb1 = i * 10
+        if s < 0.0001:
+            lb2 = i * 10
+            break
+
+    s = u"L = %d\nfyk = %d\nfck = %d\nDs = %d\nDc = %d\nLom = %d\n" % (L, fyk, fck, Ds, Dc, int(LENGTH_OF_OMOIOTHECY * 1000))
+    if isinstance(t420, PullOutExternalTassios):
+        s += "Lext = %d\n" % lext
+    s+= "-----------\n"
+    s+= "P    ds     n=1    n=3        s<0.001      s<0.0001\n"
+    s+= "420  %d : % .3f % .3f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (Ds, t420.db.s_tot[1][0], t420.db.s_tot[3][0], lb1, lb2, t420.db.s_tot[3][0] / t420.db.s_tot[1][0])
+
+    print s
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+
+    plot.plot_general(fpath, fname + "_Steel", "Length (mm)", "Stress (MPa)",
+                      (0, L/10, 0, 450),
+                      (t420.db.Ss_tot[1], "Ss = 420"),
+                      (t300.db.Ss_tot[1], "Ss = 300"),
+                      (t200.db.Ss_tot[1], "Ss = 200"),
+                      (t100.db.Ss_tot[1], "Ss = 100"))
+
+    plot.plot_general(fpath, fname + "_Concrete", "Length (mm)", "Stress (MPa)",
+                      (0, L/10, -2.5, 2.5),
+                      (t420.db.Sc_tot[1], "Ss = 420"),
+                      (t300.db.Sc_tot[1], "Ss = 300"),
+                      (t200.db.Sc_tot[1], "Ss = 200"),
+                      (t100.db.Sc_tot[1], "Ss = 100"))
+
+    plot.plot_general(fpath, fname + "_Bond", "Length (mm)", "Stress (MPa)",
+                      (0, L/10, 0, 10),
+                      (t420.db.tred_tot[1], "Ss = 420"),
+                      (t300.db.tred_tot[1], "Ss = 300"),
+                      (t200.db.tred_tot[1], "Ss = 200"),
+                      (t100.db.tred_tot[1], "Ss = 100"))
+
+    plot.plot_general(fpath, fname + "_Slip", "Length (mm)", "slip (mm)",
+                      (0, L/10, 0, 0.5),
+                      (t420.db.s_tot[1], "Ss = 420"),
+                      (t300.db.s_tot[1], "Ss = 300"),
+                      (t200.db.s_tot[1], "Ss = 200"),
+                      (t100.db.s_tot[1], "Ss = 100"))
+
+def tassios_length():
+    fck = 30
+    fyk = 400
+    P = 400
+    #L = 1000
+    Dc = 120
+    Ds = 20
+    #lext = 400
+
+    cwd = os.path.abspath(os.curdir)
+    fname = "Specimen_Length"
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+    L600 = PullOutClassicTassios(fck = fck, fyk = fyk, L = 600, Ds = Ds, Dc = Dc)
+    L600.run_recycling_positive(P, 1)
+
+    L800 = PullOutClassicTassios(fck = fck, fyk = fyk, L = 800, Ds = Ds, Dc = Dc)
+    L800.run_recycling_positive(P, 1)
+
+    L1000 = PullOutClassicTassios(fck = fck, fyk = fyk, L = 1000, Ds = Ds, Dc = Dc)
+    L1000.run_recycling_positive(P, 1)
+
+    latex = u""
+    s = u""
+    s += "fck = %d MPa\n" % fck
+    s += "fyk = %d MPa\n" % fyk
+    s += "P = %d MPa\n" % P
+    #s += "L = %d mm\n" % L
+    s += "Dc = %d mm\n" % Dc
+    s += "Ds = %d mm\n" % Ds
+    s += "Dx = %d mm\n" % int(INFINITESIMAL_LENGTH * 1000)
+    s += "Lom = %d mm \n" % int(LENGTH_OF_OMOIOTHECY * 1000)
+    s+= "-----------\n"
+    s+= "fck       n=1         n=3         s<0.001      s<0.0001\n"
+
+    for inst in (L600, L800, L1000):
+        lb1 = 0
+        lb2 = 0
+        for i, slip in enumerate(inst.db.s_tot[1]):
+            if slip < 0.001 and lb1 == 0:
+                lb1 = int(i * INFINITESIMAL_LENGTH * 1000)
+            if slip < 0.0001:
+                lb2 = int(i * INFINITESIMAL_LENGTH * 1000)
+                break
+        s1 = inst.db.s_tot[1][0]
+        s3 = inst.db.s_tot[3][0]
+        s += "   %4d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (int(inst.L * 1000),
+                                                                                    s1, s3, lb1, lb2,
+                                                                                    s3 / s1)
+        latex += "    %4d & %.7f & %.7f & %3d & %3d & %.2f\\\\\n" % (int(inst.L * 1000),
+                                                                     s1, s3, lb1, lb2,
+                                                                     s3 / s1)
+    print s
+    print latex
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+    with codecs.open(os.path.join(fpath, fname + "_latex.txt"), "w", "utf-8") as f:
+        f.write(latex)
+
+    #plot_general(fpath, fname, "Length (mm)", "slip (mm)",
+                 #(0, inst.np, 0, 0.5),
+                 #(L600.db.s_tot[1], "L=600"),
+                 #(L800.db.s_tot[1], "L=800"),
+                 #(L1000.db.s_tot[1], "L=1000"), show = True)
+
+def tassios_concrete():
+    #fck = 30
+    fyk = 500
+    P = 500
+    L = 1000
+    Dc = 120
+    Ds = 20
+    #lext = 400
+
+    cwd = os.path.abspath(os.curdir)
+    fname = "Concrete_strength"
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+    C16 = PullOutClassicTassios(fck = 16, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    C16.run_recycling_positive(500, 1)
+    C20 = PullOutClassicTassios(fck = 20, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    C20.run_recycling_positive(500, 1)
+    C25 = PullOutClassicTassios(fck = 25, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    C25.run_recycling_positive(500, 1)
+    C30 = PullOutClassicTassios(fck = 30, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    C30.run_recycling_positive(500, 1)
+    C40 = PullOutClassicTassios(fck = 40, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    C40.run_recycling_positive(500, 1)
+    C50 = PullOutClassicTassios(fck = 50, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    C50.run_recycling_positive(500, 1)
+
+    latex = u""
+    s = u""
+    #s += "fck = %d MPa\n" % fck
+    s += "fyk = %d MPa\n" % fyk
+    s += "P = %d MPa\n" % P
+    s += "L = %d mm\n" % L
+    s += "Dc = %d mm\n" % Dc
+    s += "Ds = %d mm\n" % Ds
+    s += "Dx = %d mm\n" % int(INFINITESIMAL_LENGTH * 1000)
+    s += "Lom = %d mm \n" % int(LENGTH_OF_OMOIOTHECY * 1000)
+    s+= "-----------\n"
+    s+= "fck       n=1         n=3         s<0.001      s<0.0001\n"
+
+    # Το μήκος αγκύρωσης το εξετάζω μόνο για τη διαρροή.
+    for inst in (C16, C20, C25, C30, C40, C50):
+        lb1 = 0
+        lb2 = 0
+        for i, slip in enumerate(inst.db.s_tot[1]):
+            if slip < 0.001 and lb1 == 0:
+                lb1 = int(i * INFINITESIMAL_LENGTH * 1000)
+            if slip < 0.0001:
+                lb2 = int(i * INFINITESIMAL_LENGTH * 1000)
+                break
+        s1 = inst.db.s_tot[1][0]
+        s3 = inst.db.s_tot[3][0]
+        s+= "%d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (inst.C.fck,
+                                                                               s1, s3, lb1, lb2,
+                                                                               s3 / s1)
+        latex += u"%d & %.3f & %.3f & %3d & %3d & %.2f\\\\\n" % (inst.C.fck,
+                                                                s1, s3, lb1, lb2,
+                                                                s3 / s1)
+
+    print s
+    print latex
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+    with codecs.open(os.path.join(fpath, fname + "_latex.txt"), "w", "utf-8") as f:
+        f.write(latex)
+
+
+
+def tassios_Dc():
+    fck = 30
+    fyk = 500
+    P = 500
+    L = 1000
+    #Dc = 120
+    Ds = 20
+    #lext = 400
+
+    cwd = os.path.abspath(os.curdir)
+    fname = "Concrete_diameter"
+    fpath = os.path.join(cwd, fname)
+
+    try:
+        os.mkdir(fpath)
+    except OSError:
+        pass
+
+
+    D3x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 3 * Ds)
+    D3x.run_recycling_positive(500, 1)
+    D4x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 4 * Ds)
+    D4x.run_recycling_positive(500, 1)
+    D5x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 5 * Ds)
+    D5x.run_recycling_positive(500, 1)
+    D6x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 6 * Ds)
+    D6x.run_recycling_positive(500, 1)
+    D7x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 7 * Ds)
+    D7x.run_recycling_positive(500, 1)
+    D8x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 8 * Ds)
+    D8x.run_recycling_positive(500, 1)
+    D9x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 9 * Ds)
+    D9x.run_recycling_positive(500, 1)
+    D10x = PullOutClassicTassios(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = 10 * Ds)
+    D10x.run_recycling_positive(500, 1)
+
+    latex = u""
+    s = u""
+    s += "fck = %d MPa\n" % fck
+    s += "fyk = %d MPa\n" % fyk
+    s += "P = %d MPa\n" % P
+    s += "L = %d mm\n" % L
+    #s += "Dc = %d mm\n" % Dc
+    s += "Ds = %d mm\n" % Ds
+    s += "Dx = %d mm\n" % int(INFINITESIMAL_LENGTH * 1000)
+    s += "Lom = %d mm \n" % int(LENGTH_OF_OMOIOTHECY * 1000)
+    s+= "-----------\n"
+    s+= "Dc       n=1         n=3         s<0.001      s<0.0001\n"
+
+    # Το μήκος αγκύρωσης το εξετάζω μόνο για τη διαρροή.
+
+    for inst in (D3x, D4x, D5x, D6x, D7x, D8x, D9x, D10x):
+        lb1 = 0
+        lb2 = 0
+        for i, slip in enumerate(inst.db.s_tot[1]):
+            if slip < 0.001 and lb1 == 0:
+                lb1 = int(i * INFINITESIMAL_LENGTH * 1000)
+            if slip < 0.0001:
+                lb2 = int(i * INFINITESIMAL_LENGTH * 1000)
+                break
+        s1 = inst.db.s_tot[1][0]
+        s3 = inst.db.s_tot[3][0]
+        s+= "%d : % .7f % .7f ---- lb = %3d ---- lb = %3d ---- k2 = %.2f\n" % (int(inst.Dc * 1000),
+                                                                               s1, s3, lb1, lb2,
+                                                                               s3 / s1)
+        latex += "%d & %.3f & %.3f & %3d & %3d & %.2f\\\\\n" % (int(inst.Dc * 1000),
+                                                                s1, s3, lb1, lb2,
+                                                                s3 / s1)
+
+    print s
+    print latex
+
+    with codecs.open(os.path.join(fpath, fname + ".txt"), "w", "utf-8") as f:
+        f.write(s)
+    with codecs.open(os.path.join(fpath, fname + "_latex.txt"), "w", "utf-8") as f:
+        f.write(latex[:-1])
+
+    plot_general(fpath, fname, "Length (mm)", "slip (mm)",
+                 (0, inst.np - 1, 0, 1),
+                 (D3x.db.s_tot[1], "Ds = 3xDs"),
+                 (D4x.db.s_tot[1], "Ds = 4xDs"),
+                 (D6x.db.s_tot[1], "Dc = 6xDs"),
+                 (D10x.db.s_tot[1], "Dc = 10xDs"), show = True)
+
+
+
+def store_to_shelf(fname):
+    #shelf = shelve.open(fname, flag="c", writeback=True)
+
+    shelf = {}
+
+    for fyk in (400, 500, 600):
+        sSs = "S%d" % fyk
+        for fck in (20, 30, 40):
+            sSc = "C%d" % fck
+            for Ds in (12, 16, 20):
+                sDs = "Ds%d" % Ds
+                for Dc in [d * Ds for d in (5, 6 ,8, 10)]:
+                    sDc = "Dc%d" % Dc
+
+                    key = sSs + sSc + sDs + sDc
+
+                    if key not in shelf:
+                        t = PullOutClassicTassios(fck = fck, fyk = fyk,
+                                                  L = 2000, Ds = Ds, Dc = Dc)
+                        try:
+                            t.run_recycling_positive(1.1 * fyk, 1)
+                            s1 = t.db.s_tot[1][0]
+                            s3 = t.db.s_tot[3][0]
+
+                            print s1, s3, s3/s1
+
+                            shelf[key] = t.db
+                        except ValueError:
+                            print key
+
+    with open(fname + ".cPickle", "w") as f:
+        cPickle.dump(shelf, f, protocol=cPickle.HIGHEST_PROTOCOL)
+
+
+def main():
+    a = PullOutClassicEligenhausen(fck = 30, fyk = 500, L = 1000, Ds = 25, Dc = 100)
+    a.run_recycling_positive(500, 4)
+    1
+    #plot_s(a.db)
+    #plot_Ss(a.db)
+    #plot_Sc(a.db)
+    plot_t(a.db)
+    plot_tred(a.db)
+    plot_ttred(a.db)
+    2
+if __name__ == "__main__":
+    store_to_shelf("Classic_Tassios_Fortran_LomFalse")
+
+    #store_to_shelf("Classic_Tassios_Modern_LomFalse.db")
+
+    #main()
+
+    #a = PullOutClassicEligenhausen(fck = fck, fyk = fyk, L = L, Ds = Ds, Dc = Dc)
+    #a.run_recycling_positive(fyk, 1)
+
+    #tassios_concrete()
+    #tassios_length()
+    #tassios_steel()
+    #tassios_Ds()
+    #tassios_Dc()
+    #tassios_omoiothecy_length_front()
+    #tassios_infinitesimal_length()
+
+
+
+    #a = PullOutClassicTassios(fck = 30, fyk = 800, L = 2000, Ds = 35, Dc = 105)
+    #a.run_recycling_positive(800, 1)
+
+    #s1 = a.db.s_tot[1][0]
+    #s3 = a.db.s_tot[3][0]
+
+    #print s1, s3, s3/s1
+
+
+    #a = PullOutClassicTassios(fck = 30, fyk = 220, L = 1000, Ds = 8, Dc = 100)
+    #a.run_recycling_positive(220, 1)
+
+    #s1 = a.db.s_tot[1][0]
+    #s3 = a.db.s_tot[3][0]
+
+    #print s1, s3, s3/s1
+
+    #a = PullOutClassicTassios(fck = 30, fyk = 500, L = 1000, Ds = 8, Dc = 100)
+    #a.run_recycling_positive(500, 1)
+
+    #s1 = a.db.s_tot[1][0]
+    #s3 = a.db.s_tot[3][0]
+
+    #print s1, s3, s3/s1
+
+    #######
+
+    #a = PullOutClassicTassios(fck = 30, fyk = 220, L = 1000, Ds = 18, Dc = 100)
+    #a.run_recycling_positive(220, 1)
+
+    #s1 = a.db.s_tot[1][0]
+    #s3 = a.db.s_tot[3][0]
+
+    #print s1, s3, s3/s1
+
+    #a = PullOutClassicTassios(fck = 30, fyk = 500, L = 1000, Ds = 18, Dc = 100)
+    #a.run_recycling_positive(500, 1)
+
+    #s1 = a.db.s_tot[1][0]
+    #s3 = a.db.s_tot[3][0]
+
+    #print s1, s3, s3/s1

File lines_setup.py

File contents unchanged.

File lines_test.py

File contents unchanged.
              "s" : "Slipping"}
 
 # Set Up the matplotlib figures.
-rcParams['lines.linewidth'] = 2
-rcParams['font.size'] = 14
+rcParams['lines.linewidth'] = 2.3
+rcParams['font.size'] = 12
 rcParams['font.family'] = 'serif'
 rcParams['mathtext.fontset'] = 'stix'
 rcParams['axes.labelsize'] = 'x-large'
 rcParams['xtick.labelsize'] = 'medium'
 rcParams['ytick.labelsize'] = 'medium'
-rcParams['legend.fontsize'] = 'x-large'
+rcParams['legend.fontsize'] = 'large'
 
 #import matplotlib.font_manager as fm
 #serif_font = fm.FontProperties(fname=r"/home/feanor/.fonts/Minion Pro/MinionPro-Regular.otf")
 
+def set_ax_lines_BW(ax):
+    """
+    Take each Line2D in the axes, ax, and convert the line style to be
+    suitable for black and white viewing.
+    """
+    MARKERSIZE = 3
+
+    COLORMAP = {
+        'black': {'marker': None, 'dash': (None,None), 'color' : 'black'}, # Continuous line
+        'b': {'marker': None, 'dash': [5,5], 'color' : 'b'},
+        'g': {'marker': None, 'dash': [5,3,1,3], 'color' : 'g'},
+        'r': {'marker': None, 'dash': [5,3,1,2,1,10], 'color' : 'm'},
+        'c': {'marker': None, 'dash': [5,2,5,2,5,10], 'color' : 'c'},
+        'm': {'marker': None, 'dash': [1,3], 'color' : 'r'}
+        #'y': {'marker': 'o', 'dash': [None,None]}   # Continuous line with markers
+        }
+
+    for line in ax.get_lines():
+        origColor = line.get_color()
+        #line.set_dashes(COLORMAP[origColor]['dash'])
+        line.set_marker(COLORMAP[origColor]['marker'])
+        line.set_markersize(MARKERSIZE)
+        line.set_color(COLORMAP[origColor]['color'])
+
 def _get_info(db, data):
     """
     Returns the necessary data, titles and labels for plotting and printing the
 
     if data == "Ss":
         data = db["Ss_tot"]
-        y_label = r"$S_s$ (MPa)"
+        y_label = r"$\sigma_s$ (MPa)"
         ttl = "Distribution of steel's stresses"
     elif data == "Sc":
         data = db["Sc_tot"]
-        y_label = r"$S_c$ (MPa)"
+        y_label = r"$\sigma_c$ (MPa)"
         ttl = "Distribution of concrete's stresses"
     elif data == "t":
         data = db["t_tot"]
         ttl = "Distribution of bond's stresses (reduced)"
     elif data == "es":
         data = db["es_tot"]
-        y_label = r"$e_s$ (\u2030)"
+        y_label = r"$\epsilon_s$ (\u2030)"
         ttl = "Distribution of steel's strains"
     elif data == "ec":
         data = db["ec_tot"]
-        y_label = r"$e_c$ (\u2030)"
+        y_label = r"$\epsilon_c$ (\u2030)"
         ttl = "Distribution of concrete's strains"
     elif data == "s":
         data = db["s_tot"]
     x_axis = (numpy.arange(len(data[0]), dtype = int) * int(inf_length * 1000))
 
 
-    plt.figure()
+    fig = plt.figure()
 
     # In case cycle != 0 then plot only for the given cycle
     # in case cycle == 0 then plot for all cycles.
             y_axis *= 1000
         plt.plot(x_axis, y_axis, label = "$N = %d$" %cycle)
 
-        # set the limits of the y-axis. This is not needed when all cycles are
-        # plotted because there are both negative and positive values.
-        if type_ in ("ec", "Sc"):
-            plt.ylim(ymin = min(y_axis))
-            max_ = max(y_axis)
-            if max_ < 0.0001:
-                plt.ylim(ymax = max_)
-        else:
-            plt.ylim(ymax = max(y_axis))
-            min_ = min(y_axis)
-            if min_ > -0.0001:
-                plt.ylim(ymin = min_)
         if title:
             plt.title(ttl + " for N = %d" % i)
     else:
         minimum = min([min(d) for d in data[1:]])
         maximum = max([max(d) for d in data[1:]])
 
-        if type_ in ("ec", "Sc"):
-            if maximum < 0.001:
-                plt.ylim(ymax = 0)
-        else:
-            if minimum > -0.001:
-                plt.ylim(ymin = 0)
-
         if title:
             plt.title(ttl)
 
+    # transform to black and white
+    for ax in fig.axes:
+        set_ax_lines_BW(ax)
+
+    # recreate the legend in order to get make it black and white
+    entries = []
+    for i in range(1, len(data)):
+        entries.append("$N = %d$" % i)
+    plt.legend(entries, loc = "best")
+
     plt.axhline(y=0, color = "black")
     plt.axvline(x=0, color = "black")
     plt.ylabel(y_label)
     plt.xlabel(x_label)
-    if type_ in ("ec", "Sc"):
-        plt.legend(loc = "center right")
-    else:
-        plt.legend(loc = "best")
+
     plt.grid(True)
 
+
+
 def _save(type_, path, pdf = False, png = False):
     path = os.path.join(path, FILENAMES[type_])
 
 def plot_k2(db):
     _create_plot(db, "s", cycle = 0)
     temp = ["N = %d" %1]
-    temp.extend(["N = %d" %(5 * i) for i in range(1, len(db["s_tot"]) - 1)])
+    temp.extend(["N = %d" %(4 * i + 1) for i in range(1, len(db["s_tot"]) - 1)])
 
+    plt.ylim(ymin = 0)
     plt.legend(temp, loc = "best")
     plt.show()
 
+def save_k2(db, path, pdf, png):
+    _create_plot(db, "s", cycle = 0)
+    temp = ["N = %d" %1]
+    temp.extend(["N = %d" %(4 * i + 1) for i in range(1, len(db["s_tot"]) - 1)])
+
+    plt.ylim(ymin = 0)
+    #plt.xlim(xmax = 1500)
+    plt.legend(temp, loc = "best")
+    _save("s", path, pdf = pdf, png = png)
+
 def plot_cycle(db, cycle):
     """
     Plots the distribution of all the datatypes for circle i.
 
 def _pre_vs_post_yield(db, type_):
     _create_plot(db, type_, cycle = 0)
-    plt.legend(["$P = %d$ MPa" % int(db["Ss_tot"][1][0]),
-                "$P = %d$ MPa" % int(db["Ss_tot"][2][0]),
-                "$P = %d$ MPa" % int(db["Ss_tot"][3][0])],
-               loc = "best")
+
+    entries = []
+    for i in range(1, len(db["s_tot"])):
+        entries.append("$\sigma_{s,0} = %d$ MPa" % int(db["Ss_tot"][i][0]))
+
+    plt.legend(entries, loc = "best")
 
 def plot_pre_vs_post_yield(db, type_):
     _pre_vs_post_yield(db, type_)
     # Construct the list that is going to be passed to plt.legend()
     entries = []
     for e in legend_entries:
-        entries.append("$P = %.1f$ MPa" % e)
+        entries.append("$\sigma_{s,0} = %.1f$ MPa" % e)
 
     for type_ in FILENAMES.keys():
         _pre_vs_post_yield(db, type_)
     # Construct the list that is going to be passed to plt.legend()
     entries = []
     for e in legend_entries:
-        entries.append("$P = %.1f$ MPa" % e)
+        entries.append("$\sigma_{s,0} = %.1f$ MPa" % e)
 
     for type_ in FILENAMES.keys():
         _pre_vs_post_yield(db, type_)

File pull_out_tests.pyx

+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : new.py
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+from copy import copy as ccopy
+from math import pi
+
+from scipy.optimize import brentq
+from lines import (line_slope, line_x, line_y, isascending, isdescending,
+                   interpolation, calc_area, line_intersection)
+from concrete import ModelCode2010
+from steel import B500C
+from bond import Tassios, Eligenhausen, Bond
+from constants import (INFINITESIMAL_LENGTH, LENGTH_OF_OMOIOTHECY,
+                       OMOIOTHECY_BACK, OMOIOTHECY_FRONT, LOAD_STABILITY,
+                       TOLERANCE, MAX_ITERATIONS)
+#from plot import *
+import numpy
+cimport numpy
+cimport cython
+
+class Data(object):
+    def __init__(self, L, Ds, Dc):
+        super(Data, self).__init__()
+
+        # geometrical data
+        self.L = L / 1000
+        self.Ds = Ds / 1000
+        self.Dc = Dc / 1000
+
+        self.Dx = INFINITESIMAL_LENGTH
+        self.np = int(self.L / self.Dx) + 1
+        self.L_om = LENGTH_OF_OMOIOTHECY / self.Dx
+
+        self.As = pi * self.Ds**2 / 4
+        self.Ac = pi * self.Dc**2 / 4
+        self.rho = self.As / self.Ac
+
+        np = self.np
+
+        # Lists holding the data for each node during each half-cycle.
+        self.Ss   = numpy.array([0] * np, dtype = numpy.float64)
+        self.Sc   = numpy.array([0] * np, dtype = numpy.float64)
+        self.es   = numpy.array([0] * np, dtype = numpy.float64)
+        self.ec   = numpy.array([0] * np, dtype = numpy.float64)
+        self.t    = numpy.array([0] * np, dtype = numpy.float64)
+        self.tred = numpy.array([0] * np, dtype = numpy.float64)
+        self.s    = numpy.array([0] * np, dtype = numpy.float64)
+        self.DSs  = numpy.array([0] * np, dtype = numpy.float64)
+
+        # Loaded end's stresses from previous cycles.
+        self.Ss_prev = 0.
+        self.Sc_prev = 0.
+        self.s_prev = 0.
+
+        ## Constitutive laws
+        # lists holding the points of the constitutive curve for bond and steel.
+        # each element is supposed to be a numpy array
+        self.strain = [0.] * np
+        self.stress = [0.] * np
+        self.slip = [0.] * np
+        self.bond = [0.] * np
+
+        ## min-max values for bond
+        self.smin = [0.] * np
+        self.tmin = [0.] * np
+        self.smax = [0.] * np
+        self.tmax = [0.] * np
+
+        # Database for holding test history
+        self.db = Bunch()
+
+        # Omoiothecy array
+        self._yield_reduction_omoiothecy()
+
+    def _yield_reduction_omoiothecy(self):
+        """ Calculates the reduction factor due to omoiothecy."""
+        np = self.np
+        L_om = self.L_om
+
+        self.Ko   = numpy.array([0] * np, dtype = numpy.float64)
+        for i in range(np):
+            self.Ko[i] = 1
+
+            if OMOIOTHECY_FRONT is True:
+                if i <= L_om:
+                    self.Ko[i] = (i / L_om)
+
+            if OMOIOTHECY_BACK is True:
+                if i >= np - L_om:
+                    self.Ko[i] = ((np - i - 1) / L_om)
+
+    @staticmethod
+    def db_store_data(db, Ss, Sc, es, ec, t, tred, s, DSs):
+        """
+        Appends the calculated data to the containers.
+        Should be called, after each (re)loading/unloading half-cycle.
+        """
+        assert type(db) is Bunch, "db must be a Bunch instance!!!"
+
+        db.Ss_tot.append(ccopy(Ss))
+        db.Sc_tot.append(ccopy(Sc))
+        db.es_tot.append(ccopy(es))
+        db.ec_tot.append(ccopy(ec))
+        db.t_tot.append(ccopy(t))
+        db.tred_tot.append(ccopy(tred))
+        db.s_tot.append(ccopy(s))
+        db.DSs_tot.append(ccopy(DSs))
+
+    @staticmethod
+    def db_initialize_data(db):
+        assert type(db) is Bunch, "db must be a Bunch instance!!!"
+
+        # Lists holding history of the loading history.
+        db.Ss_tot   = []
+        db.Sc_tot   = []
+        db.es_tot   = []
+        db.ec_tot   = []
+        db.t_tot    = []
+        db.tred_tot = []
+        db.s_tot    = []
+        db.DSs_tot  = []
+
+class PullOutClassic(Data):
+    def __init__(self, fck, fyk, L, Ds, Dc):
+        """
+        C : class instance (concrete)
+        S : class instance (steel)
+        B : class instance (bond)
+        """
+        super(PullOutClassic, self).__init__(L, Ds, Dc)
+
+        # Create Material Instances. Bond is a generic one.
+        self.C = ModelCode2010(fck)
+        self.S = B500C(fyk)
+        self.B = Bond()
+
+    def define_monotonic_laws(self, sign):
+        """
+        sign : int
+        must be either 1 or -1. It gives sign to the monotonic constitutive laws
+        When you start with a pull out sign == 1
+        When you start with a push in sign == -1
+        """
+        assert abs(sign) == 1, "sign must be either '1' or '-1'!!!"
+
+        self.N = 0
+        strain = sign * self.S.strain_mono
+        stress = sign * self.S.stress_mono
+        slip = sign * self.B.slip_mono
+        bond = sign * self.B.bond_mono
+
+        for i in range(self.np):
+            self.strain[i] = strain
+            self.stress[i] = stress
+            self.slip[i] = slip
+            self.bond[i] = bond
+
+    def _bound_pull(self, sf, DP):
+        """
+        Set boundary conditions for classic pull-out test.
+
+        Only the first element of each array is changed
+
+        Parameters
+        ----------
+        sf : float
+            The slip at the loaded end (mm).
+        """
+        i = 0
+
+        self.Ss[i] = self.Ss_prev + DP
+        self.Sc[i] = self.Sc_prev - DP * self.rho
+        self.es[i] = self.S.get_strain(self.strain[i], self.stress[i], self.Ss[i])
+        self.ec[i] = self.C.get_strain(self.Sc[i])
+        self.s[i]  = sf
+        self.t[i] = self.B.get_bond(self.slip[i], self.bond[i], self.s[i])
+        Ky = self.S.yield_reduction_CEB(self.es[i])
+        self.tred[i] = self.t[i] * Ky * self.Ko[i]
+        self.DSs[i] = 4 * self.tred[i] * self.Dx / self.Ds
+
+    def _bound_push(self, sf, DP):
+        """
+        Set boundary conditions for classic pull-out test.
+
+        Only the first element of each array is changed
+
+        Parameters
+        ----------
+        sf : float
+            The slip at the loaded end (mm).
+        """
+        i = 0
+
+        self.Ss[i] = self.Ss_prev + DP
+        # TODO is the boundary condition for concrete correct???
+        self.Sc[0] = self.Sc_prev
+        self.es[i] = self.S.get_strain(self.strain[i], self.stress[i], self.Ss[i])
+        self.ec[i] = self.C.get_strain(self.Sc[i])
+        self.s[i]  = sf
+        self.t[i] = self.B.get_bond(self.slip[i], self.bond[i], self.s[i])
+        Ky = self.S.yield_reduction_CEB(self.es[i])
+        self.tred[i] = self.t[i] * Ky * self.Ko[i]
+        self.DSs[i] = 4 * self.tred[i] * self.Dx / self.Ds
+
+
+    #@cython.boundscheck(False)
+    #@cython.wraparound(False)
+    #@cython.cdivision(True)
+    def _main_loop(self):
+        """
+        For a given boundary condition at the loaded end,
+        calculates the specimen, node by node.
+        """
+        cdef int i
+        cdef int np = self.np
+
+        cdef numpy.ndarray[double, ndim=1] Ss = self.Ss
+        cdef numpy.ndarray[double, ndim=1] Sc = self.Sc
+        cdef numpy.ndarray[double, ndim=1] es = self.es
+        cdef numpy.ndarray[double, ndim=1] ec = self.ec
+        cdef numpy.ndarray[double, ndim=1] s = self.s
+        cdef numpy.ndarray[double, ndim=1] t = self.t
+        cdef numpy.ndarray[double, ndim=1] tred = self.tred
+        cdef numpy.ndarray[double, ndim=1] Ko = self.Ko
+        cdef numpy.ndarray[double, ndim=1] DSs = self.DSs
+
+        cdef double rho = self.rho
+        cdef double Dx = self.Dx
+        cdef double Ds = self.Ds
+
+        cdef list strain = self.strain
+        cdef list stress = self.stress
+        cdef list slip = self.slip
+        cdef list bond = self.bond
+
+        cdef double ey = self.S.ey
+        cdef double eu = self.S.eu
+
+        temp = Dx * 500
+        red = 4 * Dx / Ds
+        for i in range(1, np):
+            Ss[i] = Ss[i-1] - DSs[i-1]
+            Sc[i] = Sc[i-1] - (Ss[i] - Ss[i-1]) * rho
+            es[i] = interpolation(stress[i], strain[i], Ss[i])
+            ec[i] = self.C.get_strain(Sc[i])
+            s[i] = (s[i-1] - ((es[i-1] + es[i]) - (ec[i-1] + ec[i])) * temp)
+            t[i] = interpolation(slip[i], bond[i], s[i])
+            #Ky = self.S.yield_reduction_CEB(self.es[i])
+            Ky = 1
+            if es[i] > ey and es[i] < eu:
+                Ky = self.S.yield_reduction_CEB(es[i])
+            tred[i] = t[i] * Ky * Ko[i]
+            DSs[i] = red * tred[i]
+
+        self.Ss = Ss
+        self.Sc = Sc
+        self.es = es
+        self.ec = ec
+        self.s = s
+        self.t = t
+        self.tred = tred
+        self.DSs = DSs
+
+    def _optimize_pull(self, sf, DP):
+        """
+        Makes the complete pull out half cycle for a given initial `sf` (slip)
+        at the loaded end.
+
+        Must be called from within `pull_load` or `pull_unload` method.
+
+        Parameters
+        ----------
+        DP : float
+            The value of stress at the loaded end (MPa).
+        sf : float
+            The slip at the loaded end (mm).
+
+        Returns
+        -------
+        self.Ss[-1] : float
+            The bar's stress at the unloaded end.
+        """
+        self._bound_pull(sf, DP)
+        self._main_loop()
+
+        return self.Ss[-1]
+
+    def _optimize_push(self, sf, DP):
+        """
+        Makes the complete push in half cycle for a given initial `sf` (slip)
+        at the loaded end.
+
+        Must be called from within `push_load` or `push_unload` method.
+
+        Parameters
+        ----------
+        sf : float
+            The slip at the loaded end (mm).
+
+        Returns
+        -------
+        self.Ss[-1] : float
+            The bar's stress at the unloaded end.
+        """
+        self._bound_push(sf, DP)
+        self._main_loop()
+
+        return self.Ss[-1]
+    def _check_convergence(self, root, output, verbose):
+        if output.converged == True:
+            if verbose == True:
+                print("--------------------------")
+                print("Convergence achieved.")
+                print("P     : {0:.1f} MPa.".format(self.Ss[0]))
+                print("slip  : {0:.10f} mm.".format(root))
+                print("Ss    : {0:.10f} MPa.".format(self.Ss[-1]))
+                print("After : {0:d} iterations".format(output.iterations))
+                print("--------------------------")
+
+            # Store boundary condition at the loaded end
+            self.Ss_prev = self.Ss[0]
+            self.Sc_prev = self.Sc[0]
+        else:
+            raise IOError("Convergence not possible.")
+
+    def pull_load(self, P, verbose = False):
+        if self.N == 0:
+            DP = P
+        else:
+            DP = P + LOAD_STABILITY
+
+        # Find the slip of balance
+        root, output = brentq(self._optimize_pull, -10, 10,
+                              args = (DP,),
+                              xtol = TOLERANCE,
+                              maxiter = MAX_ITERATIONS,
+                              full_output = True)
+
+        self._check_convergence(root, output, verbose)
+
+    def pull_unload(self, P, verbose = False):
+        DP = -P + LOAD_STABILITY
+
+        # Find the slip of balance
+        root, output = brentq(self._optimize_pull, -10, 10,
+                              args = (DP,),
+                              xtol = TOLERANCE,
+                              maxiter = MAX_ITERATIONS,
+                              full_output = True)
+
+        self._check_convergence(root, output, verbose)
+
+    def push_load(self, P, verbose = False):
+        if self.N == 0:
+            DP = -P
+        else:
+            DP = -P - LOAD_STABILITY
+
+        # Find the slip of balance
+        root, output = brentq(self._optimize_push, -10, 10,
+                              args = (DP,),
+                              xtol = TOLERANCE,
+                              maxiter = MAX_ITERATIONS,
+                              full_output = True)
+
+        self._check_convergence(root, output, verbose)
+
+    def push_unload(self, P, verbose = False):
+        DP = P - LOAD_STABILITY
+
+        # Find the slip of balance
+        root, output = brentq(self._optimize_push, -10, 10,
+                              args = (DP,),
+                              xtol = TOLERANCE,
+                              maxiter = MAX_ITERATIONS,
+                              full_output = True)
+
+        self._check_convergence(root, output, verbose)
+
+    def store_laws_and_data(self, db):
+        """ Must be overriden in the derived classes"""
+        raise NotImplementedError
+
+    @staticmethod
+    def print_summary(db):
+        for i in range(len(db.s_tot)):
+            print("N = %d \t\t slip = % .2f (mm)" %(i, db.s_tot[i][0]))
+
+        print()
+        print("k2 = %.2f" % abs(db.s_tot[3][0] / db.s_tot[1][0]))
+
+
+class PullOutClassicTassios(PullOutClassic):
+    def __init__(self, fck, fyk, L, Ds, Dc):
+        super(PullOutClassicTassios, self).__init__(fck, fyk, L, Ds, Dc)
+
+        # Overwrite Bond instance in order to make it more specific.
+        self.B = Tassios(self.C.fck, self.C.fctm)
+
+        self.db_initialize_data(self.db)
+        self.B.initialize_db(self.db)
+        self.S.initialize_db(self.db)
+
+    def store_laws_and_data(self, db):
+        """
+        Appends the new consitutive laws to the containers.
+        must be called, after each change of sign.
+        """
+        assert type(db) is Bunch, "db must be a Bunch instance!!!"
+
+        self.B.store_to_db(db, self.slip, self.bond,
+                                self.smin, self.smax)
+
+        self.S.store_to_db(db, self.strain, self.stress)
+
+        self.db_store_data(db, self.Ss, self.Sc, self.es, self.ec,
+                           self.t, self.tred, self.s, self.DSs)
+
+    def define_unloading_laws(self):
+        self.N += 1
+
+        for i in range(self.np):
+            self.smax[i] = max(self.smax[i], self.s[i])
+
+            self.slip[i], self.bond[i] = self.B.define_bond_unloading(
+                self.s[i], self.t[i], self.smin[i], self.smax[i], self.N)
+
+            self.strain[i], self.stress[i] = self.S.get_unloading_law(
+                self.es[i], self.Ss[i])
+
+    def define_reloading_laws(self):
+        self.N += 1
+
+        for i in range(self.np):
+            if i == 118:
+                pass
+
+            self.smin[i] = min(self.smin[i], self.s[i])
+
+            self.slip[i], self.bond[i] = self.B.define_bond_reloading(
+                    self.s[i], self.t[i], self.smin[i], self.smax[i], self.N)
+
+            self.strain[i], self.stress[i] = self.S.get_reloading_law(
+                self.es[i], self.Ss[i])
+
+    def run_recycling_positive(self, P, times):
+
+        self.define_monotonic_laws(1)
+        self.pull_load(0)
+        self.store_laws_and_data(self.db)
+
+        for i in range(times):
+            self.pull_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_unloading_laws()
+            self.pull_unload(P)
+            self.push_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_reloading_laws()
+            self.push_unload(P)
+
+        # finish with a pull-out
+        self.pull_load(P)
+        self.store_laws_and_data(self.db)
+
+        #self.print_summary(self.db)
+        1
+
+    def run_recycling_negative(self, P, times):
+
+        self.define_monotonic_laws(-1)
+        self.push_load(0)
+        self.store_laws_and_data(self.db)
+
+        for i in range(times):
+            self.push_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_reloading_laws()
+            self.push_unload(P)
+            self.pull_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_unloading_laws()
+            self.pull_unload(P)
+
+        # finish with a push-in
+        self.push_load(P)
+        self.store_laws_and_data(self.db)
+
+        self.print_summary(self.db)
+        1
+    def run_repeating_positive(self, P, times):
+
+        self.define_monotonic_laws(1)
+        self.pull_load(0)
+        self.store_laws_and_data(self.db)
+
+        for i in range(times):
+            self.pull_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_unloading_laws()
+            self.pull_unload(P)
+            self.store_laws_and_data(self.db)
+            self.define_reloading_laws()
+
+        # finish with a pull-out
+        self.pull_load(P)
+        self.store_laws_and_data(self.db)
+
+        self.print_summary(self.db)
+        1
+
+    def run_recycling_negative(self, P, times):
+
+        self.define_monotonic_laws(-1)
+        self.push_load(0)
+        self.store_laws_and_data(self.db)
+
+        for i in range(times):
+            self.push_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_reloading_laws()
+            self.push_unload(P)
+            self.store_laws_and_data(self.db)
+            self.define_unloading_laws()
+
+        # finish with a push-in
+        self.push_load(P)
+        self.store_laws_and_data(self.db)
+
+        self.print_summary(self.db)
+        1
+
+class PullOutExternalTassios(PullOutClassicTassios):
+    def __init__(self, fck, fyk, L, Ds, Dc, lext):
+        super(PullOutExternalTassios, self).__init__(fck, fyk, L, Ds, Dc)
+
+        assert L >= lext, "L must be greater than lext!!!"
+
+        self.lext = lext / 1000
+        self.tex = [0.] * self.np
+
+        #self.C.get_concrete_strain = self.get_concrete_strain
+
+    def _calc_external_stresses(self, P, lext):
+
+        te0 = (3/8) * (self.Ds**2 / self.Dc) * (self.S.fy / lext) * (P / self.S.fy)
+
+        lp = int(self.lext / self.Dx) + 1
+
+        for i in range(self.np):
+            if i < lp:
+                self.tex[i] = te0 * (1 - (i / lp)**2)
+
+    def _bound_pull(self, sf, DP):
+        """
+        Set boundary conditions for classic pull-out test.
+
+        Only the first element of each array is changed
+
+        Parameters
+        ----------
+        sf : float
+            The slip at the loaded end (mm).
+        """
+        i = 0
+
+        self.Ss[i] = self.Ss_prev + DP
+        self.Sc[i] = 0
+        self.es[i] = self.S.get_strain(self.strain[i], self.stress[i], self.Ss[i])
+        self.ec[i] = self.C.get_strain(self.Sc[i])
+        self.s[i]  = sf
+        self.t[i] = self.B.get_bond(self.slip[i], self.bond[i], self.s[i])
+        Ky = self.S.yield_reduction_CEB(self.es[i])
+        self.tred[i] = self.t[i] * Ky * self.Ko[i]
+        self.DSs[i] = 4 * self.tred[i] * self.Dx / self.Ds
+
+    def _bound_push(self, sf, DP):
+        """
+        Set boundary conditions for classic pull-out test.
+
+        Only the first element of each array is changed
+
+        Parameters
+        ----------
+        sf : float
+            The slip at the loaded end (mm).
+        """
+        i = 0
+
+        self.Ss[i] = self.Ss_prev + DP
+        # TODO is the boundary condition for concrete correct???
+        self.Sc[0] = 0
+        self.es[i] = self.S.get_strain(self.strain[i], self.stress[i], self.Ss[i])
+        self.ec[i] = self.C.get_strain(self.Sc[i])
+        self.s[i]  = sf
+        self.t[i] = self.B.get_bond(self.slip[i], self.bond[i], self.s[i])
+        Ky = self.S.yield_reduction_CEB(self.es[i])
+        self.tred[i] = self.t[i] * Ky * self.Ko[i]
+        self.DSs[i] = 4 * self.tred[i] * self.Dx / self.Ds
+
+    def _main_loop(self):
+        """
+        For a given boundary condition at the loaded end,
+        calculates the specimen, node by node.
+        """
+        np = self.np
+        Ss = self.Ss
+        Sc = self.Sc
+        es = self.es
+        ec = self.ec
+        s = self.s
+        t = self.t
+        tred = self.tred
+        DSs = self.DSs
+        Ko = self.Ko
+        tex = self.tex
+
+        rho = self.rho
+        Dx = self.Dx
+        Ds = self.Ds
+        Dc = self.Dc
+        temp1 = Dx * 1000 / 2
+        temp2 = 4 * Dx / Ds
+        temp3 = 4 * Dx / Dc
+
+        stress = self.stress
+        strain = self.strain
+        slip = self.slip
+        bond = self.bond
+
+        ey = self.S.ey
+        eu = self.S.eu
+
+        for i in xrange(1, np):
+            Ss[i] = Ss[i-1] - DSs[i-1]
+            self.Sc[i] = (Sc[i-1] - (Ss[i] - Ss[i-1]) * rho -
+                          tex[i] * temp3)
+            es[i] = interpolation(stress[i], strain[i], Ss[i])
+            ec[i] = self.C.get_strain(Sc[i])
+            s[i] = s[i-1] - ((es[i-1] + es[i]) - (ec[i-1] + ec[i])) * temp1
+            self.t[i] = self.B.get_bond(self.slip[i], self.bond[i], self.s[i])
+            Ky = 1
+            if es[i] > ey and es[i] < eu:
+                Ky = self.S.yield_reduction_CEB(es[i])
+            #self.tred[i] = self.t[i] * Ky * self.Ko[i]
+            tred[i] = t[i] * Ky * Ko[i]
+            DSs[i] = tred[i] * temp2
+
+        self.Ss = Ss
+        self.Sc = Sc
+        self.es = es
+        self.ec = ec
+        self.t = t
+        self.tred = tred
+        self.s = s
+        self.DSs = DSs
+
+        #for i in xrange(1, self.np):
+            #self.Ss[i] = self.Ss[i-1] - self.DSs[i-1]
+            #self.Sc[i] = (self.Sc[i-1] - (self.Ss[i] - self.Ss[i-1]) * self.rho -
+                          #self.tex[i] * 4 * self.Dx / self.Dc)
+            #self.es[i] = self.S.get_strain(self.strain[i], self.stress[i],
+                                                     #self.Ss[i])
+            #self.ec[i] = self.C.get_strain(self.Sc[i])
+            #self.s[i] = (self.s[i-1] - ((self.es[i-1] + self.es[i]) -
+                                        #(self.ec[i-1] + self.ec[i])) *
+                                        #self.Dx * 1000 / 2)
+            #self.t[i] = self.B.get_bond(self.slip[i], self.bond[i], self.s[i])
+            #Ky = self.S.yield_reduction_CEB(self.es[i])
+            #self.tred[i] = self.t[i] * Ky * self.Ko[i]
+            #self.DSs[i] = 4 * self.tred[i] * self.Dx / self.Ds
+
+    def run(self, P):
+
+        self.define_monotonic_laws(1)
+        self._calc_external_stresses(0, self.lext)
+        self.pull_load(0)
+        self.store_laws_and_data(self.db)
+
+        self._calc_external_stresses(P, self.lext)
+        self.pull_load(P)
+        self.store_laws_and_data(self.db)
+
+        1
+
+    def run_recycling_positive(self, P, times):
+
+        self.define_monotonic_laws(1)
+        self._calc_external_stresses(0, self.lext)
+        self.pull_load(0)
+        self.store_laws_and_data(self.db)
+
+        for i in range(times):
+            self._calc_external_stresses(P, self.lext)
+            self.pull_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_unloading_laws()
+            self._calc_external_stresses(0, self.lext)
+            self.pull_unload(P)
+            self._calc_external_stresses(-P, self.lext)
+            self.push_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_reloading_laws()
+            self._calc_external_stresses(0, self.lext)
+            self.push_unload(P)
+
+        # finish with a pull-out
+        self._calc_external_stresses(P, self.lext)
+        self.pull_load(P)
+        self.store_laws_and_data(self.db)
+
+        #self.print_summary(self.db)
+        1
+
+class PullOutClassicEligenhausen(PullOutClassic):
+    def __init__(self, fck, fyk, L, Ds, Dc):
+        super(PullOutClassicEligenhausen, self).__init__(fck, fyk, L, Ds, Dc)
+
+        # Overwrite Bond instance in order to make it more specific.
+        self.B = Eligenhausen(self.C.fck, self.C.fctm)
+
+        np = self.np
+
+        self.E1 = [0.] * np
+        self.E2 = [0.] * np
+        self.Esum = [0.] * np
+        self.d = [0.] * np
+        self.g = [0.] * np
+        self.Ef = [0.] * np
+        self.Efsum = [0.] * np
+        self.df = [0.] * np
+        self.tf = [0.] * np
+
+        self.db_initialize_data(self.db)
+        self.B.initialize_db(self.db)
+        self.S.initialize_db(self.db)
+
+    def store_laws_and_data(self, db):
+        """
+        Appends the new consitutive laws to the containers.
+        must be called, after each change of sign.
+        """
+        assert type(db) is Bunch, "db must be a Bunch instance!!!"
+
+        self.B.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)
+
+        self.S.store_to_db(db, self.strain, self.stress)
+
+        self.db_store_data(db, self.Ss, self.Sc, self.es, self.ec,
+                           self.t, self.tred, self.s, self.DSs)
+
+
+    def define_unloading_laws(self):
+        self.N += 1
+
+        for i in range(self.np):
+            self.smax[i], self.tmax[i], self.E1[i], self.E2[i], self.Esum[i],\
+            self.d[i], self.g[i], self.df[i], self.tf[i],\
+            self.Ef[i], self.Efsum[i], self.slip[i], self.bond[i] = \
+                self.B.define_bond_unloading(
+                    self.s[i], self.t[i], self.slip[i], self.bond[i],
+                    self.smin[i], self.smax[i], self.tmin[i],
+                    self.Esum[i], self.Efsum[i])
+
+            self.strain[i], self.stress[i] =\
+                self.S.get_unloading_law(self.es[i], self.Ss[i])
+
+    def define_reloading_laws(self):
+        self.N += 1
+
+        for i in range(self.np):
+            self.smin[i], self.tmin[i], self.E1[i], self.E2[i], self.Esum[i],\
+            self.d[i], self.g[i], self.df[i], self.tf[i],\
+            self.Ef[i], self.Efsum[i], self.slip[i], self.bond[i] = \
+                self.B.define_bond_reloading(
+                    self.s[i], self.t[i], self.slip[i], self.bond[i],
+                    self.smin[i], self.smax[i], self.tmax[i],
+                    self.Esum[i], self.Efsum[i])
+
+            self.strain[i], self.stress[i] =(
+                self.S.get_reloading_law(self.es[i], self.Ss[i]))
+
+    def run_recycling_positive(self, P, times):
+        self.define_monotonic_laws(1)
+        self.pull_load(0)
+        self.store_laws_and_data(self.db)
+
+        for i in range(times):
+            self.pull_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_unloading_laws()
+            self.pull_unload(P)
+            self.store_laws_and_data(self.db)
+            self.push_load(P)
+            self.store_laws_and_data(self.db)
+            self.define_reloading_laws()
+            self.push_unload(P)
+            self.store_laws_and_data(self.db)
+
+        # finish with a pull-out
+        self.pull_load(P)
+        self.store_laws_and_data(self.db)
+
+        self.print_summary(self.db)
+        1
+
+

File setup_pulloutests.py

+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
+
+
+#setup(
+    #cmdclass = {"build_ext": build_ext},
+    #include_dirs = [numpy.get_include()],
+    #ext_modules = [Extension("concrete",
+                             #["concrete.py"],
+                             #extra_compile_args=["-O3"]),
+                   #],
+#)
+
+#setup(
+    #cmdclass = {"build_ext": build_ext},
+    #include_dirs = [numpy.get_include()],
+    #ext_modules = [Extension("lines",
+                             #["lines.py"],
+                             #extra_compile_args=["-O3"]),
+                   #],
+#)
+
+#setup(
+    #cmdclass = {"build_ext": build_ext},
+    #include_dirs = [numpy.get_include()],
+    #ext_modules = [Extension("steel",
+                             #["steel.py"],
+                             #extra_compile_args=["-O3"]),
+                   #],
+#)
+
+#setup(
+    #cmdclass = {"build_ext": build_ext},
+    #include_dirs = [numpy.get_include()],
+    #ext_modules = [Extension("bond",
+                             #["bond.py"],
+                             #extra_compile_args=["-O3"]),
+                   #],
+#)
+
+setup(
+    cmdclass = {"build_ext": build_ext},
+    include_dirs = [numpy.get_include()],
+    ext_modules = [Extension("plot",
+                             ["plot.py"],
+                             extra_compile_args=["-O3"]),
+                   ],
+)
+
+setup(
+    cmdclass = {"build_ext": build_ext},
+    include_dirs = [numpy.get_include()],
+    ext_modules = [Extension("pull_out_tests",
+                             ["pull_out_tests.py"],
+                             extra_compile_args=["-O3"]),
+                   ],
+)

File steel_test.py

 
         show()
 
-    def run_tests(self, step, times, divisor = 1.5):
+    def run_tests(self, step, times, divisor = 1):
         """
         step : float
                It is the step of strain at each sign reversal.
+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+# module : testing.py
+##------------------------------------------------------------------------------
+
+##------------------------------------------------------------------------------
+# Author : P. Mavrogiorgos (pmav99 >at< gmail >dot< com)
+# Licence : GPL v3
+##------------------------------------------------------------------------------
+
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import os
+import shutil
+import codecs
+
+from pull_out_tests import PullOutClassicTassios
+from pull_out_tests import PullOutClassicEligenhausen
+from pull_out_tests import PullOutExternalTassios
+from pull_out_tests import PullOutClassicEligenhausenWithTassiosMonotonic
+
+import myplot
+
+def write_table_data(s_tot, path, fname ="data"):
+    lb_text = ""
+    s0_text = ""
+    k2_text = ""
+    for i in range(1, len(s_tot)):
+        lb1 = 0
+        for s in s_tot[i]:
+            if s < 0.001:
+                lb1 = s
+                break
+        lb_text += str(s_tot[i].index(lb1)) + " & "
+        s0_text += "%.2f" % s_tot[i][0] + " & "
+        k2_text += "1.00" if i == 1 else "%.2f" % (s_tot[i][0] / s_tot[1][0])
+        k2_text += " & "
+
+    lb_text = lb_text[:-2] + r"\\"
+    s0_text = s0_text[:-2] + r"\\"
+    k2_text = k2_text[:-2] + r"\\"
+
+    path = os.path.join(path, fname + ".txt")
+    with codecs.open(path, "w", "utf-8") as f:
+        f.write("\n".join([lb_text, s0_text, k2_text]))
+
+def write_latex_table(s_tot, path, fname = "table"):
+    """
+    s_tot = list
+    """
+    preamble  = r"\begin{tabular}{" + "\n"
+    preamble += r"    S[table-format = 2]" + "\n"
+    preamble += r"    S[table-format = 4]" + "\n"
+    preamble += r"    S[table-format = 1.2]" + "\n"
+    preamble += r"    S[table-format = 1.2]" + "\n"
+    preamble += r"    }" + "\n"
+
+    unit = r"\si{mm}" if len(s_tot[0]) > 1000 else r"\si{cm}"
+    temp = r"{0:10s} & {1:10s} & {2:10s} & {3:10s}\\"
+    header  = r"\toprule" + "\n"
+    header += temp.format("{$N$}",
+                          r"{$\ell_{b,1}$}",
+                          r"{$s_0$}",
+                          r"{$k_2$}") + "\n"
+    header += temp.format("{-}", unit, unit, "{-}") + "\n"
+    header += r"\midrule" + "\n"
+
+    temp = r"{0:10d} & {1:10.0f} & {2:10.2f} & {3:10s}\\"
+    body = ""
+    for i in range(1, len(s_tot)):
+        lb1 = 0
+        for s in s_tot[i]:
+            if s < 0.001: