~gcg/constrain/trunk

# -*- coding: utf-8 -*-
"""
Created on Thu Sep  3 23:17:27 2020

@author: Georg

plot stress vs. strain with Bridgeman correction

"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Button
import os
from scipy import interpolate
from scipy.optimize import curve_fit
import bottleneck as bn
import triaxiality
import configparser

cfg = configparser.ConfigParser()
cfg.read("constrain.ini")
notched_specimen = cfg.getboolean("DEFAULT", "notched_specimen")
factor = cfg.getint("DEFAULT", "nsmooth") # moving average filter -- window width for strain data
D0 = cfg.getfloat("DEFAULT", "D0") # initial diameter in cylindric gauge section
case = cfg.get("DEFAULT", "case")
dt_strain = cfg.getfloat("DEFAULT", "dt_images") # strain acquisition time period
offset_strain = cfg.getint("DEFAULT", "offset_images")  # number of images which were skipped after trigger point
nueps = cfg.getfloat("DEFAULT", "nueps") # plastic Poisson ratio

def find_limits(t1, t2):
    """ return common time range between to time axes
    """
    
    low = max(t1[0], t2[0])
    high = min(t1[-1], t2[-1])
    
    return low, high

def power_law_plasticity(strain):
    """ return power law plastic strain according to LS-Dyna material model Nr. 18"
    """
    E = 110.32
    pr = 0.33
    K = 0.6043
    n = 0.107014
    eps_ela = (E/K)**(1/(n-1))
    print("elastic strain at yield is: ", eps_ela)
    sigma = K * (eps_ela + strain)**n
    return eps_ela + strain, sigma * 1000

def load_force_LSDyna():
    """ read force from LS-Dyna axis-symmetric analysis
    """
    force_filename = "force.csv"
    data = np.genfromtxt(force_filename, delimiter=",", skip_header=2)
    time_force = data[:,0]
    force  = data[:,1] * 2 * np.pi * 1000
    dt_force = time_force[1] - time_force[0]    
    return time_force, force, dt_force

def load_force_Zwick():
    """ read force from Zwick UTM
    """
    force_filename = "force.TRA"
    data = np.genfromtxt(force_filename, delimiter=";", skip_header=16)
    time_force = data[:,0]
    force  = data[:,1]
    dt_force = time_force[1] - time_force[0]    
    return time_force, force, dt_force

def load_force_SHTB():
    """ read force from SHTB data analysis
    """
    data = np.genfromtxt("time_force_nonshifted.txt", skip_header=0)
    time_force = data[:,0]
    force  = data[:,3] * 1000 # convert kN to N
    dt_force = time_force[1] - time_force[0]
    
    #plt.plot(time_force, force)
    #plt.show()
    return time_force, force, dt_force

def load_conStrain():
    """ load strain and triaxiality data from conStrain analysis """
    
    strain_data = np.genfromtxt("strain.txt")
    indices_strain = strain_data[:,0]
    dias = strain_data[:,2]
    notch_radii = strain_data[:,3]
    
    
    eps_eq = (1./nueps) * np.log(dias[0] / dias)
    eps_eq = bn.move_mean(eps_eq,window=factor,min_count=1)
    dias = bn.move_mean(dias,window=factor,min_count=1)
    
    if notched_specimen == False:
        notch_radii = np.where(eps_eq < 0.25, 100*dias, notch_radii)
    
    notch_radii = bn.move_mean(notch_radii,window=factor,min_count=1)
    b = triaxiality.bridgman_function(dias, notch_radii)
    triaxs = triaxiality.triax_function(dias, notch_radii)
    
    
    return indices_strain, eps_eq, dias, b, triaxs

def shift_conStrain(indices_strain, nFrames, dt_strain):
    """ shift conStrain data by a number of frames to make it coincident in time with the force trigger time """
    time_strain = (indices_strain + nFrames) * dt_strain
    
    if case == "SHTB":
        time_strain += 120.0/5090 # account for strain gauge being 120 mm positioned after center of specimen
    return time_strain

def interpolate_to_common_time(time_strain):
    # construct new time axis
    dt = min(dt_force, dt_strain)
    low, high = find_limits(time_force, time_strain)
    print("common time limits: %g -- %g" % (low, high))
    n = int((high-low)/dt)
    time = np.linspace(low, high, num=n, endpoint=True)
    print("new time axis:", time[0], time[-1])
    
    # interpolate strain-based quantities onto new time axis
    f = interpolate.interp1d(time_strain, eps_eq0)
    eps_eq_i = f(time)
    f = interpolate.interp1d(time_strain, b0)
    b_i = f(time)
    f = interpolate.interp1d(time_strain, triaxs0)
    triax_i = f(time)
    f = interpolate.interp1d(time_strain, dias0)
    diameter_i = f(time)
    
    # interpolate force quantities onto new time axis
    f = interpolate.interp1d(time_force, force0)
    force_i = f(time)
    
    return eps_eq_i, diameter_i, force_i, b_i, triax_i, time



def compute_stress_measures(force, diameter, b):
    # nominal (engineering) stress
    area0 = 0.25 * np.pi * D0**2
    stress_nominal = force/area0
    
    # stress measures based on current area
    cd_mm = D0 * diameter / diameter[0]# current diameter in mm
    area = 0.25 * np.pi * cd_mm**2
    stress_cauchy = force / area
    stress_eq = stress_cauchy * b
    
    return stress_nominal, stress_cauchy, stress_eq
    

# load force file
if case == "Zwick":
    time_force, force0, dt_force = load_force_Zwick()
elif case == "SHTB":
    time_force, force0, dt_force = load_force_SHTB()
elif case == "LS-Dyna":
    time_force, force0, dt_force = load_force_LSDyna()
else:
    print("unknown force case, exiting")
    
    
print("force time axis covers %g to %g" % (time_force[0], time_force[-1]))
indices_strain, eps_eq0, dias0, b0, triaxs0 = load_conStrain()
time_strain =  shift_conStrain(indices_strain, offset_strain, dt_strain)
print("strain time axis covers %g to %g" % (time_strain[0], time_strain[-1]))  
eps_eq, diameter, force, b, triax, time = interpolate_to_common_time(time_strain)

stress_nominal, stress_cauchy, stress_eq = compute_stress_measures(force, diameter, b)

# first plot window: nominal stress and strain over time
fig, ax = plt.subplots(3, figsize=(8,7))
ax[0].plot(time, stress_nominal, color="blue", label="stress")
ax[0].legend()
ax[0].set_ylabel("nominal stress [MPa]", color="blue")
ax[0].set_xlabel("time [ms]")
ax[0].set_ylim((0,None))
ax2 = ax[0].twinx()
ax2.set_ylabel("true strain [-]", color="red") 
ax2.plot(time, eps_eq, "r-", label="strain")
ax2.set_ylim((0,None))
ax[0].grid()
ax[0].legend()


# third plot: triaxiality over time
ax[1].plot(time, triax, label="triaxiality", color="blue")
ax[1].legend()
ax[1].set_ylabel("triaxiality [-]", color="blue")
ax[1].set_xlabel("time")
ax[1].set_ylim((0, 1.2))
ax32 = ax[1].twinx()
ax32.set_ylabel("bridgman factor [-]", color="red")
ax32.plot(time, b, "r-", label=b)
ax32.set_ylim((0,1.1))
ax[1].grid()

# 3rd plot: stress / strain
l_stress_nominal, = ax[2].plot(eps_eq, stress_nominal, lw=2, label="nominal stress")
l_stress_cauchy, = ax[2].plot(eps_eq, stress_cauchy, lw=2, label="Cauchy stress")
l_stress_eq, = ax[2].plot(eps_eq, stress_eq, "g--", lw=2, label="true eq. stress, Bridgman corrected")

if case == "LS-Dyna":
    x, y = power_law_plasticity(eps_eq)
    ax[2].plot(x, y, "b--", lw=2, label="material model stress")

ax[2].grid()
ax[2].set_xlabel("equivalent strain [-]")
ax[2].set_ylabel("stress [MPa]")
ax[2].legend()
ax[2].set_xlim(0, None)
ax[2].set_ylim(0, None)





class Index(object):
    ind = 0
    def __init__(self):
        self.update_stresses()
    
    def update_stresses(self):
        offset = offset_strain + self.ind
        time_strain =  shift_conStrain(indices_strain, offset, dt_strain)
        self.eps_eq, self.diameter, self.force, self.b, self.triax, time = interpolate_to_common_time(time_strain)
        self.stress_nominal, self.stress_cauchy, self.stress_eq = compute_stress_measures(self.force, self.diameter, self.b)
        
        l_stress_nominal.set_xdata(self.eps_eq)
        l_stress_nominal.set_ydata(self.stress_nominal)
        l_stress_cauchy.set_xdata(self.eps_eq)
        l_stress_cauchy.set_ydata(self.stress_cauchy)
        l_stress_eq.set_xdata(self.eps_eq)
        l_stress_eq.set_ydata(self.stress_eq)
        ax[2].set_title("frame offset: %d, @failure: triax=%3.2f, stress=%4d, strain=%3.2f" % (offset, self.triax[-1], self.stress_eq[-1], self.eps_eq[-1]))
        plt.draw()

    def next(self, event):
        self.ind += 1
        self.update_stresses()

    def prev(self, event):
        self.ind -= 1
        self.update_stresses()

callback = Index()
axprev = plt.axes([0.2, 0.01, 0.2, 0.04])
axnext = plt.axes([0.6, 0.01, 0.2, 0.04])
bnext = Button(axnext, '+ 1 frame')
bnext.on_clicked(callback.next)
bprev = Button(axprev, '- 1 frame')
bprev.on_clicked(callback.prev)

plt.subplots_adjust(bottom=0.15, top=1)
fig.tight_layout()
plt.show()

data = np.column_stack((callback.eps_eq, callback.stress_nominal, callback.stress_cauchy, callback.stress_eq, callback.triax, callback.b))
np.savetxt("stress_strain.dat", data, header="strain, nominal_stress, Cauchy_stress, Cauchy_Bridgman_stress, triaxiality, bridgman_factor")

import sys
sys.exit()

#-----------------------
#fit Hockett-Sherby plasticity model to Cauchy-Bridgman stress
sig0 = 1500
sig1 = 500
m = 0.80
n = 2.0

bounds = [(500, 25000),
          (500, 2500),
          (0.1, 2),
          (0.1, 2.0)]
bounds = np.asarray(bounds)
bounds = bounds.T

p0 = (sig0, sig1, m, n)
indices = np.asarray(np.where(eps_eq > 0.015)).ravel()
print("indices:", indices)

if len(indices > 1):
    idx = indices[0]
else:
    idx = 0
print("start index: ", idx)
print("eps_eq:", eps_eq)

eps_pl = eps_eq[idx:]

popt, pcov = curve_fit(hocket_sherby, eps_pl, stress_eq[idx:], p0=p0, bounds=bounds)
fit = hocket_sherby(eps_pl, *popt)
format = "%10s %6s %6s %6s %6s"
format_num = "%10s %06.2f %06.2f %06.2f %06.2f"
print(format % ("parameters", "sig0", "sig1", "m", "n"))
print(format_num % ("start", *p0))
print(format_num % ("final", *popt))

plt.plot(eps_pl, stress_eq[idx:])
plt.plot(eps_pl, fit)
plt.show()


# # fit Swift-Voce plasticity model to Cauchy-Bridgman stress
# A = 1501
# eps0 = 0.001
# n = 0.40
# k0 = 668
# Q = 1093
# b = 1.54
# alpha = 0.5

# bounds = [(1000, 2000),
#           (0.001, 0.2),
#           (0.1,0.5),
#           (200,1000),
#           (200,2000),
#           (0.1, 100),
#           (0,1)]
# bounds = np.asarray(bounds)
# bounds = bounds.T

# p0 = (A, eps0, n, k0, Q, b, alpha)
# idx = np.where(eps_eq > 0.0025)[0][0]
# eps_pl = eps_eq[idx:]
# print("start index: ", idx)
# popt, pcov = curve_fit(swift_voce, eps_pl, stress_eq[idx:], p0=p0, bounds=bounds)
# fit = swift_voce(eps_pl, *popt)
# format = "%10s %6s %6s %6s %6s %6s %6s %6s"
# format_num = "%10s %06.2f %06.2f %06.2f %06.2f %06.2f %06.2f %06.2f"
# print(format % ("parameters", "A", "eps0", "n", "k0", "Q", "b", "alpha"))
# print(format_num % ("start", *p0))
#print(format_num % ("final", *popt))

# --------------