Commit a811c580 authored by Michael Rudolf's avatar Michael Rudolf
Browse files

Resolve #2

 - Removed redundant scripts
 - Changed conversion scripts to use `rstevaluation` for conversion.
 - Changed how filenames are generated for plots.
parent 7e4630a2
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 13 15:13:56 2018
@author: Warsitzka, Michael
This script can be used to convert tdms files of the velocity stepping test into text files and to plot the vst-data.
"""
#%%=============IMPORT=====================================================
import os
import numpy as np
from nptdms import TdmsFile
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from os.path import splitext
import itertools
import csv
#%%==============PARAMETERS=================================================
A = 0.022619 # area of shear zone (= surface area of lid) (m^2)
li = 0.0776 # inner lever (center cell to center of shear zone)
lo = 0.1250 # outer lever (center cell to hinge point)
le= 1.611655493 # load point level
v= 30 # shear velocity (mm/min)
#g = 9.81 # acceleration due to gravity
normalstress= 2000 # normal stress in Pa, only for .txt or asc. files
#%%==================NAMES=================================================
path_in = 'input/' #define path of input data
path_out = 'output/' #define path of input data
curve_fit= 'logx' #define type of curve fit for shear velocity-shear stress curve:
#>>'lin' = linear:
# [slope, intercept] = logfit(x,y,'linear')
# yApprox = (intercept)+(slope)*x
#>> 'logx' = logarithmic:
# [slope, intercept] = logfit(x,y,'logx')
# yApprox = (intercept)+(slope)*log10(x)
#>>'logy' = exponential:
# [slope, intercept] = logfit(x,y,'logy')
# yApprox = (10^intercept)*(10^slope).^x
#>>'loglog' = powerlaw:
# [slope, intercept] = logfit(x,y,'loglog')
# yApprox = (10^intercept)*x.^(slope)
#%%==========================READ & CONVERT DATA==============================================
files = [i for i in os.listdir(path_in) \
if i.endswith(".tdms") and \
('VST' in i or 'vst' in i)]
n = len(files)
names = [splitext(files[i])[0] for i in range(0,n)]
dfvst_raw = np.empty(n, dtype=object)
dfvst = np.empty(n, dtype=object)
i=0
for i in range(0,n):
tdms = TdmsFile(path_in + files[i])
groups = tdms.groups() # get group names
channels = tdms.group_channels(groups[0]) #get channel names
inc = list((tdms.object(groups[0], 'Velocity').properties).values())[2] # time increment between each point
dfvst_raw[i] = tdms.object(groups[0]).as_dataframe() #load tdms as data frame
shearforce = dfvst_raw[i]['Shear Force'] #velocity in [N]
normalforce = dfvst_raw[i]['Normal Force'] #velocity in [N]
liddispl = dfvst_raw[i]['Lid Displacement'] #velocity in [mm]
vel = dfvst_raw[i]['Velocity'] # shear velocity in [mm/s]
if filter(lambda x: x < 0, vel) == True:
vel[vel < 0] = nan
time = np.arange(0, len(vel)*inc, inc) # time in [s]
displ = (np.cumsum(inc*vel)) #cumulative displacement in [mm]
sigma = normalforce/A #calculate normal stress
tau = (shearforce*lo)/(li*A) #calculate shear stress
fric = (tau/sigma) #friction []
dil = -(liddispl-liddispl[0]) # calculation dilation
dfvst[i]= pd.DataFrame.from_dict({'time': time,
'displacement': displ,
'velocity': vel,
'normalstress': sigma,
'shearstress': tau,
'friction': fric,
'dilation':dil}, orient='index').transpose()
#%%==========FIT SHEAR LOAD CURVE============================================
if curve_fit== 'lin': #'linear'
logvel=vel
logfric=fric
elif curve_fit== 'logx': #'logarithmic'
logvel=np.log10(vel)
logfric=fric
elif curve_fit== 'logy': #'exponential'
logfric=np.log10(fric)
logvel=vel
elif curve_fit== 'loglog': #'powerlaw'
logvel=np.log10(vel)
logfric=np.log10(fric)
# linear fit and evaluation:
slope, intercept,r,p,std = stats.linregress(logvel,logfric)
fitfric= np.polyval((slope,intercept),logvel)
# Standard deviations:
m=len(logvel)
s = np.sum((fric-intercept-slope*logvel)**2)/(m-2)
std_slope = np.sqrt(m*s/(m*np.sum(logvel**2)-np.sum(logvel)**2)) #standard error Coef. friction
std_intercept = np.sqrt(np.sum(logvel**2) * s/ (m*np.sum(logvel**2)-np.sum(logvel)**2)) #standard error cohesion
# rescale the approximation results for plotting
if curve_fit == 'lin':
name_fit='y = '+str(round(slope,3))+'*x +'+str(round(intercept,3))
elif curve_fit=='logx':
name_fit = r'y = ({:.4f}$\pm${:.4f})log10(x) + ({:.4f}$\pm${:.4f})'.format(slope,
std_slope,intercept,std_intercept)
if curve_fit=='logy':
fitfric=10**fitfric
logfric=10**logfric
name_fit = 'y = (10^('+str(round(intercept,3))+')*10^('+str(round(slope,3))+')^x'
elif curve_fit=='loglog':
ran=10**ran
fitfric=10**fitfric
logfric=10**logfric
name_fit= 'y = 10^('+str(round(intercept,3))+')*x^('+str(round(slope,3))+')'
# Calculate MSE and R2
MSE = np.mean((logfric-fitfric)**2) # mean squared error.
Covy = np.cov(fitfric, fric) # = cov(estimated values, experimental values)
R2 = (Covy[1]**2)/(np.var(fitfric)*np.var(fric)) #variance
#%%=======================SAVE VST DATA TO TXT FILE=================================================
# data_vst = pd.DataFrame({'Time [s]': time, 'Shear velocity [mm/s]': vel, 'Normal force [N]': normalforce, 'Shear force [N]': shearforce})
data_vst = np.array((time, vel, normalforce, shearforce))
#data_vst.to_csv(path+names[i] + '.txt',index=None, sep='\t', mode='a')
header_vst = 'Time [s]'+'\t' 'Shear velocity [mm s^-1]'+ '\t'+'Normal force [N]'+'\t'+'Shear force [N]'
with open(path_out + names[i] + '.txt', 'w') as f:
f.write(header_vst + '\n')
for row in data_vst.T:
csv.writer(f, delimiter='\t', lineterminator='\n').writerow(row)
f.close
# -*- coding: utf-8 -*-
"""
Converts all *.asc files in the folder given in `path_in` and creates a time
series plot in the folder `path_out`. This script requires `rstevaluation` to
be installed on the system.
Created on Tue Sep 25 14:54:43 2018
@author: analab
@author: M.Warsitzka, M.Rudolf
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import csv
import pylab
from scipy import stats
import os
import fnmatch
import glob, shutil
import pickle
from os.path import isfile, join, basename, splitext
from nptdms import TdmsFile
import collections
import itertools
import codecs
#%%==========================PARAMETERS & VARIABLES==============================================
A = 0.022619 # area of shear zone (= surface area of lid) (m^2)
li = 0.0776 # inner lever (center cell to center of shear zone)
lo = 0.1250 # outer lever (center cell to hinge point)
le= 1.611655493 #load point level
v= 30 #shear velocity (mm/min)
g = 9.81 # acceleration due to gravity
prec = 125 # rounding precision for normal load in file names
velnoise = 1e-2
stressnoise = 0.025
smoothwindow = 51 # define window for smooth shear stress curve
projectname = '324-01_QuartzneusandBern_basal'
TSname= projectname+'_ts'
path_in = 'input/'
path_out = 'output/'
#%%==========================READ & CONVERT DATA==============================================
for file in os.listdir(path_in):
files,n = glob.glob1(path_in,"*.asc"),len(glob.glob1(path_in,'*.asc'))
DF = np.empty(n, dtype=object)
len_m = np.zeros(n)
Sigma_avg = np.zeros(m)
for i in range(0,len(files)):
if files[i].endswith('.asc'): # loop over all asc files
name = splitext(files[i])[0] # extract names of asc files
with codecs.open(path_in + files[i], encoding='utf-8-sig') as f:
data_load = np.loadtxt(f) # load file for further calculations
# write data to new varables with meaningfull names
time = data_load[:,0]
shearload = data_load[:,1] #in kg
velocity = data_load[:,2]
dilation = data_load[:,3]
normalload = data_load[:,4] #in kg
# Corrections
if time[0] >= 0:
time = time - time[0] # set starting time to 0
else:
time = time
dilation = -(dilation-dilation[0]) # set initial dilation to 0
shearforce = shearload *g #in N
normalforce = normalload *g #in N
sigma = normalforce/A # kg to normal stress in [Pa]
tau = (shearforce*lo)/(li*A) # kg to shear stress in [Pa]
velocity = velocity/60 # transform to [mm/s]
velocity[velocity < 0] = 0
inc = time[1]-time[0] # timesteps, first one = 0
#incdisp = np.diff(time)*velocity[:-1] # incremental displacement in mm
displacement = np.cumsum(inc*velocity[:-1]) # cummulative displacement in mm
Sigma_avg[i] = np.average(sigma[sigma>0]) # only for the file name
DF[i]= pd.DataFrame.from_dict({'time': time,
'displacement': displacement,
'velocity': velocity,
'normalforce': normalforce,
'normalstress': sigma,
'shearforce': shearforce,
'shearstress': tau,
'dilation':dilation}, orient='index').transpose()
# open(path+ str(name)+'_'+str(int(avgsigma))+'.txt', 'w').close() #write and save calculated data in new txt file
# with open(path+ str(name)+'_'+str(int(avgsigma))+'.txt', 'r+') as output_set:
# write_data = csv.writer(output_set, delimiter='\t')
# write_data.writerow(['Time', 'Shear stress', 'Velocity', 'Dilation','Normal stress', 'Displacement'])
# write_data.writerows(zip(time,tau,velocity,dilation,sigma,displacement))
# output_set.closed
# current_header= np.genfromtxt((path_in+str(files_txt[i])),
# delimiter='\t',
# dtype='unicode')[0]
len_m[i] = len(time)
print(name +' processed')
break
TS = np.zeros((n+1,maxlen_df))#np.zeros((n+1, maxlen_df)) #array of zeros for time series of all files
#%%==========================EVALUATE STRESSES==============================================
for i in range(0,len(DF)):
current_set = DF[i]
#create df of time series for txt file:
current_ts = np.concatenate(([pd.Series(current_set.time).values],
[pd.Series(current_set.shearforce).values]))
# ts_dspl = np.zeros(ts.shape)
# for l in range(0, current_ts.shape[0]):
# ts_dspl[l,0] = np.copy(ts[l,0] * vel/60)
# ts_dspl[l,1] = np.copy(ts[l,1])
# l = l+1
#convert load (kg) to stress (Pa)
# for k in range(1,current_ts.shape[1]):
# for l in range(0,current_ts.shape[0]):
# ts_dspl[l,k] = ts[l,k]*9.81*loadlev/A
# l = l+1
# k= k+1
current_ts = np.asarray([np.pad(a, (0, maxlen_df - len(a)), 'constant', constant_values=np.NaN) for a in current_ts])
TS[i+1] = current_ts[1]
if i == len_m_max_ind:
TS[0] = current_ts[0]
#%%=======================PLOT TS=================================================
plt.rcParams['figure.figsize'] = (10, 8)
linecolor= ['','red', 'orange','gold', 'green', 'royalblue', 'k']
t = int(m/3)
fig3=plt.figure()
for i in range(0,t):
sigma_avg = int((data_peak[:,0][i]+data_peak[:,0][i+t]++data_peak[:,0][i+2*t])/3)
plt.plot(TS.T[:,0], np.zeros(len(TS.T[:,0])), linewidth= 0.5, color=linecolor[i+1], label=str(int(sigma_avg))+' Pa')
plt.plot(TS.T[:,0], TS.T[:,i*3+1:(i+1)*3+1], linewidth= 0.5, color=linecolor[i+1])
plt.legend(fontsize=8,
facecolor='w',
edgecolor='k',
framealpha=1,
loc= 'upper right',
title="Normal stress")
plt.xlabel('Shear displacement [mm]')
plt.ylabel(r'Shear stress $\tau$ [Pa]')
plt.xlim(0, max(TS.T[:,0]))
fig3.suptitle(rstname, y=0.92)
plt.savefig(path_out + TSname, bbox_inches='tight', edgecolor='w')
plt.close()
#%% Save ts data to txt files:
sigma_str = "\t".join([str(int(x.round(-1))) for x in Sigma_avg])#.split(',')
header_ts = 'Time [s]'+'\t'+'Normal stress [Pa]:'+'\t'+sigma_str
with open(path_out + TSname + '.txt', 'w') as f:
f.write(header_ts + '\n')
for row in TS.T:
csv.writer(f, delimiter='\t', lineterminator='\n').writerow(row)
f.close
\ No newline at end of file
import numpy as np
import rstevaluation.tools.files as rstfiles
import rstevaluation.tools.plots as rstplots
def main():
projectname = '393-01_01_Foamglass100'
path_in = 'Legacy/FileConversion/input'
path_out = 'Legacy/FileConversion/output'
# =========== PARAMETERS & VARIABLES ============
var = {
'A': 0.022619, # area of shear zone (= surface area of lid) (m^2)
'li': 0.0776, # inner lever (center cell to center of shear zone)
'lo': 0.1250, # outer lever (center cell to hinge point)
'le': 1.611655493, # load point level
'v': 30, # shear velocity (mm/min)
'g': 9.81, # acceleration due to gravity
'prec': 125, # rounding precision for normal load in file names
'velnoise': 1e-2,
'stressnoise': 0.025,
'smoothwindow': 51, # define window for smooth shear stress curve
'nsamples': 10**6 # Maximum number of samples for downsampling
}
exp_data = [
rstfiles.convert(path_in, f, var)
for f in os.listdir(path_in)
if f.endswith('.asc')
]
normal_stresses = [
'%i' % (np.floor(np.mean(d['normalstress'])/var['prec'])*var['prec'])
for d in exp_data
]
rstplots.plotts(path_out, projectname, exp_data, normal_stresses)
rstfiles.saveTS(path_out, projectname, exp_data, normal_stresses)
if __name__ == '__main__':
main()
# -*- coding: utf-8 -*-
"""
Converts all *.tdms files in the folder given in `path_in` and creates a time
series plot in the folder `path_out`. This script requires `rstevaluation` to
be installed on the system.
Created on Tue Sep 25 14:54:43 2018
@author: analab
@author: M.Warsitzka, M.Rudolf
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import csv
import pylab
from scipy import stats
import os
import fnmatch
import glob, shutil
import pickle
from os.path import isfile, join, basename, splitext
from nptdms import TdmsFile
import collections
import itertools
#%%==========================PARAMETERS & VARIABLES==============================================
A = 0.022619 # area of shear zone (= surface area of lid) (m^2)
li = 0.0776 # inner lever (center cell to center of shear zone)
lo = 0.1250 # outer lever (center cell to hinge point)
le= 1.611655493 #load point level
v= 30 #shear velocity (mm/min)
g = 9.81 # acceleration due to gravity
prec = 125 # rounding precision for normal load in file names
velnoise = 1e-2
stressnoise = 0.025
smoothwindow = 51 # define window for smooth shear stress curve
params_conv = [A,li,lo,g,prec]
params_eval = [velnoise, stressnoise, A, v, le]
projectname = '393-01_01_Foamglass100'
TSname= projectname+'_ts'
path_in = 'input/'
path_out = 'output/'
#%%==========================READ & CONVERT DATA==============================================
for file in os.listdir(path_in):
files,n = glob.glob1(path_in,"*.tdms"), len(glob.glob1(path_in,'*.tdms'))
DF = np.empty(n, dtype=object)
len_m = np.zeros(n)
for i in range(0,n):
name = splitext(files[i])[0]
tdms = TdmsFile(path + files[i])
groups = tdms.groups()
channels = tdms.group_channels(groups[0])
# time increment between each point:
inc = list((tdms.object(groups[0], 'Velocity').properties).values())[2]
df_raw=tdms.object(groups[0]).as_dataframe()
Length[i] = df_raw.shape[0]
shearforce = df_raw['Shear Force'] #Shear force in [N]
normalforce = df_raw['Normal Force'] #normal forcein [N]
liddispl = df_raw['Lid Displacement'] #velocity in [mm]
velocity = df_raw['Velocity'] #velocity in [mm/s]
#convert colums:
time = np.arange(0, Length[i]*inc, inc) # time in [s]
displacement = np.cumsum(inc*velocity) #cumulative displacement in [mm]
sigma = normalforce/A
tau = (shearforce*lo)/(li*A)
dilation = -(liddispl-liddispl[0])
df[i]= pd.DataFrame.from_dict({'time': time,
'displacement': displacement,
'velocity': velocity,
'normalforce': normalforce,
'normalstress': sigma,
'shearforce': shearforce,
'shearstress': tau,
'dilation':dilation}, orient='index').transpose()
print(name +' processed')
len_m[i] = len(time)
break
len_m_max = int(max(len_m))
len_m_max_ind = next(i for i, j in enumerate(len_m) if j == max(len_m))
TS = np.zeros((n+1,maxlen_df))#np.zeros((n+1, maxlen_df)) #array of zeros for time series of all files
#%%==========================EVALUATE STRESSES==============================================
Sigma_avg = np.zeros(m)
for i in range(0,len(DF)):
current_set = DF[i]
#create df of time series for txt file:
current_ts = np.concatenate(([pd.Series(current_set.time).values],
[pd.Series(current_set.shearforce).values]))
# ts_dspl = np.zeros(ts.shape)
# for l in range(0, current_ts.shape[0]):
# ts_dspl[l,0] = np.copy(ts[l,0] * vel/60)
# ts_dspl[l,1] = np.copy(ts[l,1])
# l = l+1
#convert load (kg) to stress (Pa)
# for k in range(1,current_ts.shape[1]):
# for l in range(0,current_ts.shape[0]):
# ts_dspl[l,k] = ts[l,k]*9.81*loadlev/A
# l = l+1
# k= k+1
sigma = pd.Series(current_set.shearstress).values
Sigma_avg[i] = np.average(sigma[sigma>0])
current_ts = np.asarray([np.pad(a, (0, maxlen_df - len(a)), 'constant', constant_values=np.NaN) for a in current_ts])
TS[i+1] = current_ts[1]
if i == len_m_max_ind:
TS[0] = current_ts[0]
#%%=======================PLOT TS=================================================
plt.rcParams['figure.figsize'] = (10, 8)
linecolor= ['','red', 'orange','gold', 'green', 'royalblue', 'k']
t = int(m/3)
fig3=plt.figure()
for i in range(0,t):
normalstress_avg = int((data_peak[:,0][i]+data_peak[:,0][i+t]++data_peak[:,0][i+2*t])/3)
plt.plot(TS.T[:,0], np.zeros(len(TS.T[:,0])), linewidth= 0.5, color=linecolor[i+1], label=str(int(normalstress_avg))+' Pa')
plt.plot(TS.T[:,0], TS.T[:,i*3+1:(i+1)*3+1], linewidth= 0.5, color=linecolor[i+1])
plt.legend(fontsize=8,
facecolor='w',
edgecolor='k',
framealpha=1,
loc= 'upper right',
title="Normal stress")
plt.xlabel('Shear displacement [mm]')
plt.ylabel(r'Shear stress $\tau$ [Pa]')
plt.xlim(0, max(TS.T[:,0]))
fig3.suptitle(rstname, y=0.92)
plt.savefig(path_out + TSname, bbox_inches='tight', edgecolor='w')
plt.close()
#%% Save ts data to txt files:
sigma_str = "\t".join([str(int(x.round(-1))) for x in Sigma_avg])#.split(',')
header_ts = 'Time [s]'+'\t'+'Normal stress [Pa]:'+'\t'+sigma_str
with open(path_out + TSname + '.txt', 'w') as f:
f.write(header_ts + '\n')
for row in TS.T:
csv.writer(f, delimiter='\t', lineterminator='\n').writerow(row)
f.close
\ No newline at end of file
import numpy as np
import rstevaluation.tools.files as rstfiles
import rstevaluation.tools.plots as rstplots
def main():
projectname = '393-01_01_Foamglass100'
path_in = 'Legacy/FileConversion/input'
path_out = 'Legacy/FileConversion/output'
# =========== PARAMETERS & VARIABLES ============
var = {
'A': 0.022619, # area of shear zone (= surface area of lid) (m^2)
'li': 0.0776, # inner lever (center cell to center of shear zone)
'lo': 0.1250, # outer lever (center cell to hinge point)
'le': 1.611655493, # load point level
'v': 30, # shear velocity (mm/min)
'g': 9.81, # acceleration due to gravity
'prec': 125, # rounding precision for normal load in file names
'velnoise': 1e-2,
'stressnoise': 0.025,
'smoothwindow': 51, # define window for smooth shear stress curve
'nsamples': 10**6 # Maximum number of samples for downsampling
}
exp_data = [
rstfiles.convert(path_in, f, var)
for f in os.listdir(path_in)
if f.endswith('.tdms')
]
normal_stresses = [
'%i' % (np.floor(np.mean(d['normalstress'])/var['prec'])*var['prec'])
for d in exp_data
]
rstplots.plotts(path_out, projectname, exp_data, normal_stresses)
rstfiles.saveTS(path_out, projectname, exp_data, normal_stresses)
if __name__ == '__main__':
main()
# -*- coding: utf-8 -*-
"""
Converts all *.tdms files in the folder given in `path_in` and creates a
velocity stepping plot for each of them in the folder `path_out`. This script
requires `rstevaluation` to be installed on the system.
@author: M.Rudolf
"""
import os
import numpy as np
import rstevaluation.tools.analysis as rstanalysis
import rstevaluation.tools.files as rstfiles
import rstevaluation.tools.plots as rstplots
def main():
projectname = '393-01_01_Foamglass100'
path_in = 'Legacy/FileConversion/input'
path_out = 'Legacy/FileConversion/output'
# =========== PARAMETERS & VARIABLES ============
var = {
'A': 0.022619, # area of shear zone (= surface area of lid) (m^2)
'li': 0.0776, # inner lever (center cell to center of shear zone)
'lo': 0.1250, # outer lever (center cell to hinge point)
'le': 1.611655493, # load point level
'v': 30, # shear velocity (mm/min)
'g': 9.81, # acceleration due to gravity
'prec': 125, # rounding precision for normal load in file names
'velnoise': 1e-2,
'stressnoise': 0.025,
'smoothwindow': 51, # define window for smooth shear stress curve
'nsamples': 10**6 # Maximum number of samples for downsampling
}
exp_data = [
rstfiles.convert(path_in, f, var)
for f in os.listdir(path_in)
if f.endswith('.tdms')
]
normal_stresses = [
'%i' % (np.floor(np.mean(d['normalstress'])/var['prec'])*var['prec'])
for d in exp_data
]
for exp in exp_data:
pfit, perr, name_fit = rstanalysis.vst_analysis(exp)
rstplots.plotVST(path_out, exp, pfit, perr, name_fit)
rstfiles.saveTS(path_out, projectname, exp_data, normal_stresses)