ANNEXE: Python script for the analytic solution
==================================================


.. code-block:: text

    ################################################

    # JOINT_MECA_ENDO: Shear test #

    ################################################

    from math import*

    Import numpy as nP

    import scipy.optimize as opt

    import matplotlib.pyplot as plt

    ###############################################################

    SELECTION DES DONNEES FROM ENTREE:

    ###############################################################

    # Fashion = ['Ilaria', 'CIH']

    fashion = 'CIH'

    if mode == 'Ilaria':

        # Friction coefficient

        mu = 1. # Excerpt from ssnap162a line 35

        Bn = 137622240.236444 # Excerpt from lcejdm line 166

        Bt = 32281760.05546218 # Excerpt from lcejdm line 167

        # Energy dissipated during damage:

        D1 = 1. # Excerpt from lcejdm line 174

    if mode == 'CIH':

        # Cohesion

        C = 1.4e6 # Email extract from CIH on 8/02

        #C = 1.e5 # Excerpt from Ilaria's test

        # Maximum shear stress

        T = 1.1e6 # Email extract from CIH on 8/02

        #T = 0.9*C # Excerpt from Ilaria's test

        # Friction coefficient

        mu = 0.95 # Email extract from CIH on 8/02

        # Normal stiffness:

        Kn = 1.e12 # Excerpt from ssnap162a line 30

        # Tangential stiffness:

        Set = 2. *Kn # Excerpt from ssnp162a line 32

        # Residual cohesion

        c_bar = 0.6

        # Constant compressive stress:

        sigma_n_comp = 0.4e6

        # Peak shape parameters

        m1 = 3. # Excerpt from lcejdm line 172

        m2 = 0.5 # Excerpt from lcejdm line 173

        # No integration (alpha from 0 to 1 per step step)

        step=0.0001

        # Max damage (between 0 and 1)

        endo_max = 0.90

    ###############################################################

    # Definition of S and R functions as a function of alpha:

    ###############################################################

    # Definition of the S function (alpha):

    def S (alpha):

        return ((1-alpha) **m1)/(alpha**m2)

    # Definition of the function derived from S (alpha):

    def Sp (alpha):

        return ((m2-m1) *alpha-m2)* ((1-alpha) **(m1-1))/(alpha** (m2+1))

    # Definition of the second derivative of S (alpha):

    def S2p (alpha):

        return ((m1-m2) *(m1-m2-1.)*alpha**2+2.*m2*(m1-m2-1.)*alpha+
                m2*(m2+1.))* (1.-alpha**(m1-2)) /alpha** (m2+2)

    # Definition of the R function (alpha):

    def R (alpha):

        return (2). *(Sp (alpha))**2)/(S2p (alpha))) -S (alpha)

    #############################################################

    # Finding the max of R (alpha):

    #############################################################

    def findMaxr ():

        val_max = opt.fmin (lambda alpha: -R (alpha), 1.e-10)

        return R (val_max [0])

    #############################################################

    # Calculating Bn, Bt, C, and T:

    #############################################################

    def calcBN ():

        Bn = (-Sp (alpha_max) *T**2)/(2. *D1*S (alpha_max) **2)

        Return Bn

    def calcBT ():

        Bt = (-Sp (alpha_max) *(C**2-mu**2*t**2))/(2.*D1*S (alpha_max)**2)

        Return Bt

    def calCC ():

        C = sqrt (2. *D1* (Bn*(mu**2) +Bt))) *S (alpha_max)/(sqrt (-Sp (alpha_max))))

        Return C

    def calculation ():

        T = sqrt (2. *D1*Bn) *S (alpha_max) /sqrt (-Sp (alpha_max))

        Return T

    #############################################################

    # Calculation of delta and sigma as a function of alpha:

    #############################################################

    def join_meca_endo (Kn, Kt, Bn, Bt, Bn, Bt, m1, m2, mu, c_bar, D1, alpha):

        delta_t = (sqrt (2. *D1/ (- (Mu**2*BN+Bt)*Sp (alpha)))) *(Kt+ (Mu**2*BN+Bt)*S (alpha)) /Kt+

        (mu*sigma_n_comp+c_bar) /Kt)

        sigma_t = sqrt (2. *D1* (mu**2*bn+bt)) *S (alpha) /sqrt (-Sp (alpha)) +mu*sigma_n_comp+c_bar

        delta_n = -sigma_n_comp/kn+mu*sqrt (2.*D1/ (- (Mu**2*BN+Bt) *Sp (alpha)))

        delta_n_trac = sqrt (2*D1/ (-Bn*Sp (alpha)))) *(Kn+Bn*S (alpha)) /Kn+C_bar/ (Mu*kn)

        sigma_n = S (alpha) *sqrt (2.*D1*Bn) /sqrt (-Sp (alpha)) +c_bar/mu

        return delta_t, sigma_t, delta_n, delta_n_trac, sigma_n

    #############################################################

    # Printing the curves:

    #############################################################

    def plotCurve (x, y, title, label_x, label_x, label_y, output_folder_path, filename):

        plt.figure (title)

        plt.plot (x, y, label=title)

        plt.xlabel (label_x)

        plt.ylabel (label_y)

        plt.suptitle (title)

        plt.legend ()

        plt.savefig (output_folder_path+filename, dpi=300)

        plt.close ()

    if __name__ == '__main__':

        # Alpha_max calculation:

        alpha_max = (-m2+sqrt (m1*m2/ (m1-m2+1))))/(m1-m2)

        # Search for as much R as possible:

        R_max = findMaxR ()

        # Parameter calculations based on input data:

        if mode == 'Ilaria':

            C = calCC ()

            T = calCT ()

        if mode == 'CIH':

            D1 = max ((max (T**2/ (2.*Kn), C**2/ (2.*Kt)) *(-Sp (alpha_max)/(S (alpha_max)**2)) *R_max*10.) ,1.) ,1.)

            Bn = calcBN ()

            Bt = CalcBt ()

            # D_min calculation:

        D_min = max (T**2/ (2.*Kn), C**2/ (2.*Kt)) *(-Sp (alpha_max)/(S (alpha_max)**2)) *R_max

        # Data printing:

        Print '\n'

        print '===== VALEUR DES PARAMETRES ===== '

        print '---------------------------------------------'

        Print 'Kn:', Kn, '\ t\ t\ t Kt:', Kt

        print 'C:', C, '\ t\ t\ t:', T

        print 'mu: ', mu,'\ t\ t\ t c_bar: ', c_bar

        print 'Bn: ', Bn,'\ t\ t Bt: ', Bt

        print 'alpha_max: ', alpha_max,'\ t R_max: ', R_max:', R_max

        print 'D1: ', D1,'\ t\ t sigma_n_comp: ', sigma_n_comp

        print 'm1: ', m1,'\ t\ t\ t m2: ', m2

        # Parameter validation tests:

        Print '\n'

        print '===== TEST BY VALIDATION DES PARAMETRES ===== '

        print '---------------------------------------------'

        print 'Test D1 > D_min:\ t\ t\ t', D1>=d_min

        If not D1>D_min:

            Print '/! \ D1 = ', D1, 'and D_min =', D_min

            print 'Test Mu*t < C < sqrt (2.)*Mu*t:\ t\ t', Mu*t < C < sqrt (2.) *Mu*t

        if not mu*T < C < sqrt (2.)*Mute*:

            Print '/! \ Mu*t = ', Mu*t, 'C =', C, 'sqrt (2.) *mu*T = ', sqrt (2.) *Mu*t

            print 'Test Kn > BN*R_Max:\ t\ t\ t', Kn>BN*R_Max

        If not KN>BN*R_Max:

            Print '/! \ Kn = ', Kn, 'BN*R_max =', BN*R_max = ', BN*R_max

            print 'Test Kt > BT*R_Max:\ t\ t\ t', kt>BT*R_Max

        If not KT>BT*R_Max:

            Print '/! \ Kt = ', Kt, 'BT*R_Max =', BT*R_Max

            Print 'Test Bn > 0\ t\ t\ t\ t', Bn>0.

        If not Bn>0. :

            Print '/! \ Bn = ', Bn

            print 'Test Bt > 0\ t\ t\ t\ t', Bt>0.

        If not Bt>0. :

            Print '/! \ Bt = ', Bt

        print '---------------------------------------------'

        Print '\n'

        # Delta and sigma calculations:

        list_delta_t = [:ref:`0. <0.>`]

        list_sigma_t = [:ref:`0. <0.>`]

        list_delta_n = [:ref:`0. <0.>`]

        list_delta_n_trac = [:ref:`0. <0.>`]

        list_sigma_n = [:ref:`0. <0.>`]

        for alpha in np.arange (0. +step, endo_max+step, step):

            delta_t, sigma_t, delta_n, delta_n_trac, sigma_n = joint_meca_endo (Kn, Kt, Bn, Bt, m1, m2, m2, m2, mu, c_bar, D1, alpha)

            list_delta_t.append (delta_t*1e3)

            list_sigma_t.append (sigma_t/1.e6)

            list_delta_n.append (delta_n*1e3)

            list_delta_n_trac.append (delta_n_trac*1e3)

            list_sigma_n.append (sigma_n/1.e6)