import csv
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib import gridspec
from pylab import genfromtxt
from scipy.signal import argrelextrema
import pandas as pd
from matplotlib.ticker import LogLocator
from pylab import genfromtxt
from matplotlib.lines import Line2D
from itertools import cycle

plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
#%%
# -----------------------------------------------------------------------------
# Color definitions (RGB)
# Color schemes generated with "colobrewer2.org"
# -----------------------------------------------------------------------------
myBlue     = ('#eff3ff','#c6dbef','#9ecae1','#6baed6','#3182bd','#08519c') ; myBlue = myBlue[::-1]
myGreen    = ('#006d2c','#31a354','#74c476','#a1d99b','#c7e9c0','#edf8e9')
myOrange   = ('#a63603','#e6550d','#fd8d3c','#fdae6b','#fdd0a2','#feedde')
myOrange2  = ('#ffffd4','#fee391','#fec44f','#fe9929','#d95f0e','#993404') ; myOrange2 = myOrange2[::-1]
myRed      = ('#fee5d9','#fcbba1','#fc9272','#fb6a4a','#de2d26','#a50f15') ; myRed = myRed[::-1]
Diverging6 = ('#8c510a','#d8b365','#f6e8c3','#c7eae5','#5ab4ac','#01665e')
Diverging8 = ('#8c510a','#bf812d','#dfc27d','#f6e8c3','#c7eae5','#80cdc1','#35978f','#01665e')
Diverging9 = ('#d73027','#f46d43','#fdae61','#fee090','#ffffbf','#e0f3f8','#abd9e9','#74add1','#4575b4')
Diverging10 = ('#a50026','#d73027','#f46d43','#fdae61','#fee090','#e0f3f8','#abd9e9','#74add1','#4575b4','#313695')
Diverging11 = ('#a50026','#d73027','#f46d43','#fdae61','#fee090','#ffffbf','#e0f3f8','#abd9e9','#74add1','#4575b4','#313695')

#%%            
# -----------------------------------------------------------------------------
# Define matplotlib(fonts) dictionaries
families = ['serif', 'sans-serif', 'cursive', 'fantasy', 'monospace']
weights = ['light', 'normal', 'medium', 'semibold', 'bold', 'heavy', 'black']
Linestyles = [':',':',':',':','-.','-.']
marker = ['v', 's', '^', 'd', '>', '<']

# Costumize the appearence
# ------------------------
# Font size and display dpi
plt.rcParams['figure.dpi'] = 200
plt.rcParams['savefig.dpi'] = 400

# Fonts
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = families[0]      # Set the font family, available are: families = ['serif', 'sans-serif', 'cursive', 'fantasy', 'monospace']
plt.rcParams['font.weight'] = weights[1]       # Set the fonts weights, available are: ['light', 'normal', 'medium', 'semibold', 'bold', 'heavy', 'black']

# Text
plt.rcParams['figure.titlesize'] = 'xx-small'
plt.rcParams['text.usetex'] = False             # Enable conversion to Latex font
plt.rcParams['text.color'] = 'black'           # Textcolor (inside the graph)
           
# Axes
plt.rcParams['axes.linewidth'] = 0.5
plt.rcParams['xtick.color'] = 'black'           # Textcolor of the x-axes
plt.rcParams['ytick.color'] = 'black'           # Textcolor of the y-axes
plt.rcParams['axes.edgecolor'] = 'black'      # Textcolor of the y-axes
plt.rcParams['axes.labelcolor'] = "black"

# Grid
plt.rcParams['grid.linewidth'] = 0.3        # Set the grid color
plt.rcParams['grid.color'] = 'silver'          # Set the grid color

# Legend
plt.rcParams['legend.fontsize'] = 'small'
#plt.rcParams['legend.framealpha'] = None       # No transparent legend backgroun
plt.rcParams['legend.framealpha'] = 0.75       
plt.rcParams['legend.fancybox'] = False
plt.rcParams['legend.facecolor'] = 'white'
#plt.rcParams['legend.handlelength'] = 1

#%%
# -----------------------------------------------------------------------------
# Functions / Definitions
# -----------------------------------------------------------------------------
def cm2inch(*tupl):
    inch = 2.54
    if isinstance(tupl[0], tuple):
        return tuple(i/inch for i in tupl[0])
    else:
        return tuple(i/inch for i in tupl)
    
def moving_average(interval, window_size):
    window= np.ones(int(window_size))/float(window_size)
    return np.convolve(interval, window, 'same')

def read_TCUI_IncDr(filepath):
    """Load vertical strain, deviatoric stress and mean effective stress from a TMU output file."""
    data = np.loadtxt(filepath, skiprows=1, usecols=(2,8,9,10))
    q = data[:,2]-data[:,1]; eps1 = -data[:,0]*100
    p = (-data[:,1]-data[:,2]-data[:,3])/3
    return q, eps1, p  # returns (q, eps1, p)

def read_TCUI_numgeo(filepath):
    """Load vertical strain, deviatoric stress and mean effective stress from a TMU output file."""
    data = np.loadtxt(filepath, skiprows=1, usecols=(1,2,3,8))
    q = data[:,0]-data[:,1]; eps1 = -data[:,3]*100
    p = (-data[:,0]-data[:,1]-data[:,2])/3
    return q, eps1, p  # returns (q, eps1, p)

#%%
# picture settings
nfigs = 2
height = 12; width = 15  # cm

#--------------------------------------------------
#%%

files_TCUI_IncDr = [('IncDr/TCUI19.out',  myOrange[1])]
files_TCUI_numgeo = [('numgeo/TCUI/TCUI19-print-out/EALL_element_1.dat',  myOrange[1])]
# files_TCUI_numgeo = [('numgeo/TCUI/TCUI19_2-print-out/EALL_element_1.dat',  myOrange[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs)
plt.subplot(gs[0])
for filepath, color in files_TCUI_IncDr:
    q, eps1, p = read_TCUI_IncDr(filepath)
    plt.plot(p, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TCUI_numgeo:
    q, eps1, p = read_TCUI_numgeo(filepath)
    plt.plot(p, q, color=color, ls='--', lw=1.25)
legend_handles = [Line2D([0], [0], color=myOrange[1], ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color=myOrange[1], ls='--', lw=2, label='numgeo, Forward-Euler') ]
leg1 = plt.legend(handles=legend_handles, loc='best')
plt.gca().add_artist(leg1)

cc_values = ["0\\%"]
cc_handles = []
for (filepath, color), cc in zip(files_TCUI_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))
plt.legend(handles=cc_handles, loc='upper left')

plt.xlim(0.,120.)
plt.xticks(np.arange(0, 120.1, 30), fontsize=14)
plt.ylim(-80.,80.)
plt.yticks(np.arange(-80., 80.1, 20), fontsize=14)
plt.xlabel(r"Mean effective stress, p$'$ in kPa", fontsize=14)
plt.ylabel('Deviatoric stress, q in kPa', fontsize=14) 
    
plt.subplot(gs[1])
for filepath, color in files_TCUI_IncDr:
    q, eps1, p = read_TCUI_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TCUI_numgeo:
    q, eps1, p = read_TCUI_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)
 
legend_handles = [Line2D([0], [0], color=myOrange[1], ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color=myOrange[1], ls='--', lw=2, label='numgeo, Forward-Euler') ]
leg1 = plt.legend(handles=legend_handles, loc='best')
plt.gca().add_artist(leg1)

cc_handles = []
for (filepath, color), cc in zip(files_TCUI_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))
plt.legend(handles=cc_handles, loc='upper left')

plt.xlim(-8.,6.)
plt.xticks(np.arange(-8, 6.1, 2), fontsize=14)
plt.ylim(-80.,80.)
plt.yticks(np.arange(-80., 80.1, 20), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, q in kPa', fontsize=14) 
     
plt.show()
# fig.savefig('TCUI19.pdf', bbox_inches='tight')
fig.savefig('TCUI19.png', bbox_inches='tight')

#--------------------------------------------------
#%%
files_TCUI_IncDr = [('IncDr/TCUI19_BC.out',  myOrange[1])]
files_TCUI_numgeo = [('numgeo/TCUI/TCUI19_BC-print-out/EALL_element_1.dat',  myOrange[1])]
# files_TCUI_numgeo = [('numgeo/TCUI/TCUI19_2-print-out/EALL_element_1.dat',  myOrange[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs)
plt.subplot(gs[0])
for filepath, color in files_TCUI_IncDr:
    q, eps1, p = read_TCUI_IncDr(filepath)
    plt.plot(p, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TCUI_numgeo:
    q, eps1, p = read_TCUI_numgeo(filepath)
    plt.plot(p, q, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color=myOrange[1], ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color=myOrange[1], ls='--', lw=2, label='numgeo, Forward-Euler') ]
leg1 = plt.legend(handles=legend_handles, loc='best')
plt.gca().add_artist(leg1)

cc_values = ["1\\%"]
cc_handles = []
for (filepath, color), cc in zip(files_TCUI_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))
plt.legend(handles=cc_handles, loc='upper left')

plt.xlim(0.,120.)
plt.xticks(np.arange(0, 120.1, 30), fontsize=14)
plt.ylim(-80.,80.)
plt.yticks(np.arange(-80., 80.1, 20), fontsize=14)
plt.xlabel(r"Mean effective stress, p$'$ in kPa", fontsize=14)
plt.ylabel('Deviatoric stress, q in kPa', fontsize=14) 
    
plt.subplot(gs[1])
for filepath, color in files_TCUI_IncDr:
    q, eps1, p = read_TCUI_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TCUI_numgeo:
    q, eps1, p = read_TCUI_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)
 
legend_handles = [Line2D([0], [0], color=myOrange[1], ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color=myOrange[1], ls='--', lw=2, label='numgeo, Forward-Euler') ]
leg1 = plt.legend(handles=legend_handles, loc='best')
plt.gca().add_artist(leg1)

cc_handles = []
for (filepath, color), cc in zip(files_TCUI_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))
plt.legend(handles=cc_handles, loc='upper left')

plt.xlim(-8.,6.)
plt.xticks(np.arange(-8, 6.1, 2), fontsize=14)
plt.ylim(-80.,80.)
plt.yticks(np.arange(-80., 80.1, 20), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, q in kPa', fontsize=14) 
     
plt.show()
# fig.savefig('TCUI19.pdf', bbox_inches='tight')
fig.savefig('TCUI19_BC.png', bbox_inches='tight')

#%%

files_TCUI_IncDr = [('IncDr/TCUI2.out',  myOrange[1])]
files_TCUI_numgeo = [('numgeo/TCUI/TCUI2-print-out/EALL_element_1.dat',  myOrange[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs)  
plt.subplot(gs[0])
for filepath, color in files_TCUI_IncDr:
    q, eps1, p = read_TCUI_IncDr(filepath)
    plt.plot(p, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TCUI_numgeo:
    q, eps1, p = read_TCUI_numgeo(filepath)
    plt.plot(p, q, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color=myOrange[1], ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color=myOrange[1], ls='--', lw=2, label='numgeo, Forward-Euler') ]
leg1 = plt.legend(handles=legend_handles, loc='best')
plt.gca().add_artist(leg1)

cc_values = ["6\\%"]
cc_handles = []
for (filepath, color), cc in zip(files_TCUI_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))
plt.legend(handles=cc_handles, loc='upper left')

plt.xlim(0.,125.)
plt.xticks(np.arange(0, 125.1, 25), fontsize=14)
plt.ylim(-150.,150.)
plt.yticks(np.arange(-150., 150.1, 50), fontsize=14)
plt.xlabel(r"Mean effective stress, p$'$ in kPa", fontsize=14)
plt.ylabel('Deviatoric stress, q in kPa', fontsize=14) 
    
plt.subplot(gs[1])
for filepath, color in files_TCUI_IncDr:
    q, eps1, p = read_TCUI_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TCUI_numgeo:
    q, eps1, p = read_TCUI_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color=myOrange[1], ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color=myOrange[1], ls='--', lw=2, label='numgeo, Forward-Euler') ]
leg1 = plt.legend(handles=legend_handles, loc='best')
plt.gca().add_artist(leg1)

cc_handles = []
for (filepath, color), cc in zip(files_TCUI_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))
plt.legend(handles=cc_handles, loc='upper left')

plt.xlim(-5.,5.)
plt.xticks(np.arange(-5, 5.1, 2.5), fontsize=14)
plt.ylim(-150.,150.)
plt.yticks(np.arange(-150., 150.1, 50), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, q in kPa', fontsize=14) 
 
plt.show()
# fig.savefig('TCUI20.pdf', bbox_inches='tight')
fig.savefig('TCUI2.png', bbox_inches='tight')

