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

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')
PurpleBlue = ('#762a83','#9970ab','#c2a5cf','#e7d4e8','#d9f0d3','#a6dba0','#5aae61','#1b7837')
#%%            
# -----------------------------------------------------------------------------
# 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'] = 500
plt.rcParams['savefig.dpi'] = 500

# 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_oed_IncDr(filepath):
    """Load vertical stress and void ratio from an OED output file."""
    data = np.loadtxt(filepath, skiprows=1, usecols=(8, 14))
    s1=-data[:, 0]; e=data[:, 1]
    return s1,e  # returns (s1, e)

def read_TMD_IncDr(filepath):
    """Load vertical strain, volumetric strain and deviatoric stress from a TMD output file."""
    data = np.loadtxt(filepath, skiprows=1, usecols=(2,3,4,8,9))
    q = data[:,4]-data[:,3]; eps1 = -data[:,0]*100
    epsv = (-data[:,0]-data[:,1]-data[:,2])*100
    return q, eps1, epsv  # returns (q, eps1, epsv)

def read_TMU_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_oed_numgeo(filepath):
    """Load vertical stress and void ratio from an OED output file."""
    data = np.loadtxt(filepath, skiprows=1, usecols=(2, 13))
    s1=-data[:, 0]; e=data[:, 1]
    return s1,e  # returns (s1, e)

def read_TMD_numgeo(filepath):
    """Load vertical strain, volumetric strain and deviatoric stress from a TMD output file."""
    data = np.loadtxt(filepath, skiprows=1, usecols=(1,2,7,8,9))
    q = data[:,0]-data[:,1]; eps1 = -data[:,3]*100
    epsv = (-data[:,2]-data[:,3]-data[:,4])*100
    return q, eps1, epsv  # returns (q, eps1, epsv)

def read_TMU_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 = 14  # cm

"""Drained Monotonic Triaxial tests""" 
files_TMD_IncDr = [
    ('IncDr/TMD6.out',  myOrange[1]),
    ('IncDr/TMD11.out',  myGreen[1]),
    ('IncDr/TMD16.out',  myBlue[1])]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD6-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD11-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD16-print-out/EALL_element_1.dat',  myBlue[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='white', ls='-', lw=2, label=r'$c_c$ = 0%'),
    Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
    Line2D([0], [0], color='black', ls='--', lw=2, label='numgeo, Forward-Euler') ]
plt.legend(handles=legend_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,1500.)
plt.yticks(np.arange(0., 1500.1, 300.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
    
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)
  
legend_handles = [Line2D([0], [0], color='white', ls='-', lw=2, label=r'$c_c$ = 0%'),
    Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
    Line2D([0], [0], color='black', ls='--', lw=2, label='numgeo, Forward-Euler') ]
plt.legend(handles=legend_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-4,2.)
plt.yticks(np.arange(-4., 2.1, 1.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Untreated1.png', bbox_inches='tight')

###### -----------------------------------------------------------------------------
files_TMD_IncDr = [
    ('IncDr/TMD5.out',  myOrange[1]),
    ('IncDr/TMD10.out',  myGreen[1]),
    ('IncDr/TMD15.out',  myBlue[1]),]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD5-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD10-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD15-print-out/EALL_element_1.dat',  myBlue[1]),]
    
fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])   
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='white', ls='-', lw=2, label=r'$c_c$ = 0%'),
    Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
    Line2D([0], [0], color='black', ls='--', lw=2, label='numgeo, Forward-Euler') ]
plt.legend(handles=legend_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,1500.)
plt.yticks(np.arange(0., 1500.1, 300.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
        
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)

for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='white', ls='-', lw=2, label=r'$c_c$ = 0%'),
    Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
    Line2D([0], [0], color='black', ls='--', lw=2, label='numgeo, Forward-Euler') ]
plt.legend(handles=legend_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-2,2.)
plt.yticks(np.arange(-2., 2.1, 1), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Untreated2.png', bbox_inches='tight')

###### -----------------------------------------------------------------------------
files_TMD_IncDr = [
    ('IncDr/TMD1.out',  myOrange[1]),
    ('IncDr/TMD6.out',  myGreen[1]),
    ('IncDr/TMD20.out',  myBlue[1]),]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD1-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD6-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD20-print-out/EALL_element_1.dat',  myBlue[1]),]
  
fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])   
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='white', ls='-', lw=2, label=r'$c_c$ = 0%'),
    Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
    Line2D([0], [0], color='black', ls='--', lw=2, label='numgeo, Forward-Euler') ]
plt.legend(handles=legend_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,500.)
plt.yticks(np.arange(0., 500.1, 100.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
        
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)

for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='white', ls='-', lw=2, label=r'$c_c$ = 0%'),
    Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
    Line2D([0], [0], color='black', ls='--', lw=2, label='numgeo, Forward-Euler') ]
plt.legend(handles=legend_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-6,2.)
plt.yticks(np.arange(-6., 2.1, 2), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Untreated3.png', bbox_inches='tight')

###### -----------------------------------------------------------------------------
files_TMD_IncDr = [
    ('IncDr/TMD6.out',  myOrange[1]),
    ('IncDr/TMD7.out',  myGreen[1]),
    ('IncDr/TMD8.out',  myBlue[1]),
    ('IncDr/TMD9.out', Diverging6[1])]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD6-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD7-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD8-print-out/EALL_element_1.dat',  myBlue[1]),
    ('numgeo/TMD/TMD9-print-out/EALL_element_1.dat', Diverging6[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])   
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color='black', ls='--', lw=2, label='numgeo, Forward-Euler') ]

cc_handles = []
cc_values = ["0\\%", "2.8\\%", "6.3\\%", "11.5\\%"]
for (filepath, color), cc in zip(files_TMD_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))

all_handles = legend_handles + cc_handles
plt.legend(handles=all_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,2000.)
plt.yticks(np.arange(0., 2000.1, 400.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
        
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)

for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color='black', 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_TMD_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='lower left')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-8,2.)
plt.yticks(np.arange(-8.,2.1, 2), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Treated1.png', bbox_inches='tight')

###### -----------------------------------------------------------------------------
files_TMD_IncDr = [
    ('IncDr/TMD11.out',  myOrange[1]),
    ('IncDr/TMD12.out',  myGreen[1]),
    ('IncDr/TMD13.out',  myBlue[1]),
    ('IncDr/TMD14.out', Diverging6[1])]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD11-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD12-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD13-print-out/EALL_element_1.dat',  myBlue[1]),
    ('numgeo/TMD/TMD14-print-out/EALL_element_1.dat', Diverging6[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])   
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

cc_handles = []
cc_values = ["0\\%", "4.6\\%", "7.2\\%", "11.1\\%"]
for (filepath, color), cc in zip(files_TMD_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))

all_handles = legend_handles + cc_handles
plt.legend(handles=all_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,2400.)
plt.yticks(np.arange(0., 2400.1, 600.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
        
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)

for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color='black', 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_TMD_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='lower left')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-6,2.)
plt.yticks(np.arange(-6., 2.1, 2), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Treated2.png', bbox_inches='tight')

###### -----------------------------------------------------------------------------
files_TMD_IncDr = [
    ('IncDr/TMD16.out',  myOrange[1]),
    ('IncDr/TMD17.out',  myGreen[1]),
    ('IncDr/TMD18.out',  myBlue[1]),
    ('IncDr/TMD19.out', Diverging6[1])]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD16-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD17-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD18-print-out/EALL_element_1.dat',  myBlue[1]),
    ('numgeo/TMD/TMD19-print-out/EALL_element_1.dat', Diverging6[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])   
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

cc_handles = []
cc_values = ["0\\%", "6.8\\%", "8\\%", "13.9\\%"]
for (filepath, color), cc in zip(files_TMD_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))

all_handles = legend_handles + cc_handles
plt.legend(handles=all_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,6000.)
plt.yticks(np.arange(0., 6000.1, 1000.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
        
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)

for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color='black', 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_TMD_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='lower left')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-5,2.)
plt.yticks(np.arange(-5., 2.1, 1), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Treated3.png', bbox_inches='tight')

###### -----------------------------------------------------------------------------
files_TMD_IncDr = [
    ('IncDr/TMD1.out',  myOrange[1]),
    ('IncDr/TMD2.out',  myGreen[1]),
    ('IncDr/TMD3.out',  myBlue[1]),
    ('IncDr/TMD4.out', Diverging6[1])]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD1-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD2-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD3-print-out/EALL_element_1.dat',  myBlue[1]),
    ('numgeo/TMD/TMD4-print-out/EALL_element_1.dat', Diverging6[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])   
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

cc_handles = []
cc_values = ["0\\%", "3.6\\%", "7.7\\%", "12.8\\%"]
for (filepath, color), cc in zip(files_TMD_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))

all_handles = legend_handles + cc_handles
plt.legend(handles=all_handles, loc='best')
plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,1200.)
plt.yticks(np.arange(0., 1200.1, 200.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
        
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)

for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)

legend_handles = [Line2D([0], [0], color='gray', ls='-', lw=2, label='IncDr, Forward-Euler'),
Line2D([0], [0], color='black', 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_TMD_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='lower left')
plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-4,2.)
plt.yticks(np.arange(-4., 2.1, 1), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Treated4.png', bbox_inches='tight')

###### -----------------------------------------------------------------------------
files_TMD_IncDr = [
    ('IncDr/TMD20.out',  myOrange[1]),
    ('IncDr/TMD21.out',  myGreen[1]),
    ('IncDr/TMD22.out',  myBlue[1]),
    ('IncDr/TMD23.out', Diverging6[1])]

files_TMD_numgeo = [
    ('numgeo/TMD/TMD20-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMD/TMD21-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMD/TMD22-print-out/EALL_element_1.dat',  myBlue[1]),
    ('numgeo/TMD/TMD23-print-out/EALL_element_1.dat', Diverging6[1])]

fig = plt.figure(figsize=cm2inch(nfigs*width,height+0.5))
gs = gridspec.GridSpec(1,nfigs, height_ratios=[1], width_ratios = [1]*nfigs, wspace=0.25)
plt.subplot(gs[0])   
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, q, color=color, ls='--', lw=1.25)

cc_handles = []
cc_values = ["0\\%", "1.4\\%", "5.1\\%", "9.2\\%"]
for (filepath, color), cc in zip(files_TMD_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))

all_handles = legend_handles + cc_handles
plt.legend(handles=all_handles, loc='best')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(0.0,2400.)
plt.yticks(np.arange(0., 2400.1, 600.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
        
plt.subplot(gs[1])
for filepath, color in files_TMD_IncDr:
    q, eps1, epsv = read_TMD_IncDr(filepath)
    plt.plot(eps1, epsv, color=color, ls='-', lw=1.25, alpha=0.6)

for filepath, color in files_TMD_numgeo:
    q, eps1, epsv = read_TMD_numgeo(filepath)
    plt.plot(eps1, epsv, color=color, ls='--', lw=1.25)

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

cc_handles = []
for (filepath, color), cc in zip(files_TMD_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='lower left')

plt.xlim(0.,30.)
plt.xticks(np.arange(0., 30.1, 5.), fontsize=14)
plt.ylim(-12,2.)
plt.yticks(np.arange(-12., 2.1, 2), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel(r'Volumetric strain, $\varepsilon_v$ in %', fontsize=14)
plt.axhline(y=0, color='k', ls=(0, (8, 6)), lw=0.5)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMD_Treated5.png', bbox_inches='tight')

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

"""UnDrained Monotonic Triaxial tests""" 
files_TMU_IncDr = [
    ('IncDr/TMU1.out',  myOrange[1]),
    ('IncDr/TMU2.out',  myGreen[1]),
    ('IncDr/TMU3.out',  myBlue[1]),
    ('IncDr/TMU4.out',  Diverging6[1]),
    ('IncDr/TMU5.out',  PurpleBlue[1]),
    ('IncDr/TMU6.out',  myRed[1])]

files_TMU_numgeo = [
    ('numgeo/TMU/TMU1-print-out/EALL_element_1.dat',  myOrange[1]),
    ('numgeo/TMU/TMU2-print-out/EALL_element_1.dat',  myGreen[1]),
    ('numgeo/TMU/TMU3-print-out/EALL_element_1.dat',  myBlue[1]),
    ('numgeo/TMU/TMU4-print-out/EALL_element_1.dat',  Diverging6[1]),
    ('numgeo/TMU/TMU5-print-out/EALL_element_1.dat',  PurpleBlue[1]),
    ('numgeo/TMU/TMU6-print-out/EALL_element_1.dat',  myRed[1])]

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

cc_handles = []
cc_values = ["0\\%", "8.4\\%", "0\\%", "7.6\\%", "0\\%", "7.7\\%"]
for (filepath, color), cc in zip(files_TMU_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))

all_handles = legend_handles + cc_handles
plt.legend(handles=all_handles, loc='best')

plt.xlim(0.,10.)
plt.xticks(np.arange(0., 10.1, 2.), fontsize=14)
plt.ylim(0.0,4000.)
plt.yticks(np.arange(0., 4000.1, 500.), fontsize=14)
plt.xlabel(r'Vertical strain, $\varepsilon_1$ in %', fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)
    
plt.subplot(gs[1])
for filepath, color in files_TMU_IncDr:
    q, eps1, p = read_TMU_IncDr(filepath)
    plt.plot(p, q, color=color, ls='-', lw=1.25, alpha=0.6)
    
for filepath, color in files_TMU_numgeo:
    q, eps1, p = read_TMU_numgeo(filepath)
    plt.plot(p, q, color=color, ls='--', lw=1.25)

cc_handles = []
cc_values = ["0\\%", "8.4\\%", "0\\%", "7.6\\%", "0\\%", "7.7\\%"]
for (filepath, color), cc in zip(files_TMU_IncDr, cc_values):
    cc_handles.append(Line2D([0], [0], color=color, lw=2, linestyle='-', label=rf'$c_c = {cc}$'))

all_handles = legend_handles + cc_handles
plt.legend(handles=all_handles, loc='best')

plt.xlim(0.,1500.)
plt.xticks(np.arange(0., 1500.1, 300.), fontsize=14)
plt.ylim(0.0,4000.)
plt.yticks(np.arange(0., 4000.1, 500.), fontsize=14)
plt.xlabel(r"Mean effective stress, p$'$ in kPa", fontsize=14)
plt.ylabel('Deviatoric stress, $q$ in kPa', fontsize=14)

plt.show()
# fig.savefig('TMD.pdf', bbox_inches='tight')
fig.savefig('TMU.png', bbox_inches='tight')
