Commit 5c5f9b2a by Michael Rudolf

### Added all scripts by @warsitzka_at_ig.cas.cz

parent 245625c9
 # -*- 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 -*- """ Created on Tue Sep 25 14:54:43 2018 @author: analab """ 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
 # -*- coding: utf-8 -*- """ Created on Tue Sep 25 14:54:43 2018 @author: analab """ 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
 # -*- coding: utf-8 -*- """ Created on Mon Oct 22 10:16:04 2018 @author: analab """ #%%===========================IMPORT========================================== 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 nptdms import TdmsFile import collections import itertools import RST_convert_old projectname = '419-01_Quartzsand_Zaragosa_L-7080S' path_in = 'input/' path_out = 'output/' name_ts= projectname+'_ts' #%%==========================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) v= 30 #shear velocity (mm/min) prec = 125 # rounding precision for normal load in file names params_conv = [A,li,lo,prec] #%%==========================READ & CONVERT DATA============================================== def convert_tdms(path,params): files,N = glob.glob1(path,"*.tdms"), len(glob.glob1(path,'*.tdms')) A,li,lo,prec = params[0],params[1],params[2],params[3] # values for calculation DF = np.empty(N, dtype=object) Length = np.zeros(N) Sigma_avg = 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_N_max = int(max(len_m)) len_N_max_ind = next(i for i, j in enumerate(len_m) if j == max(len_m)) DFTS = np.zeros((N+1,len_N_max))#np.zeros((n+1, maxlen_df)) #array of zeros for time series of all files for i in range(0,len(DF)): current_set = DF[i] current_ts = np.concatenate(([pd.Series(current_set.time).values], [pd.Series(current_set.shearforce).values])) current_ts = np.asarray([np.pad(a, (0, len_N_max- len(a)), 'constant', constant_values=np.NaN) for a in current_ts]) DFTS[i+1] = current_ts[1] if i == len_N_max_ind: DFTS[0] = current_ts[0] return DFTS #%%=======================PLOT DFTS================================================= def RST_plot(out_path, TSname, DFTS): # plt.rcParams['savefig.format'] = 'pdf' # plt.rcParams['figure.figsize'] = (14, 10) # plt.rcParams['font.size'] = 10 # plt.rcParams['font.family'] = 'Arial Unicode MS' # plt.rcParams['savefig.dpi'] = 1000 # plt.rcParams['figure.figsize'] = (10, 8) # for k in range(0,DFTS.shape[0]): # plt.plot(DFTS.T[:,0], DFTS.T[:,k], linewidth= 0.5, color='k') # k= k+1 # plt.xlabel('Shear displacement [mm]') # plt.ylabel(r'Shear stress [Pa] ') # plt.xlim(0, max(DFTS[0])) # plt.ylim(0, 40) # plt.savefig(out_path + TSname, bbox_inches='tight', edgecolor='w') # plt.close() Sigma = [500,1000,2000,4000,8000,16000,500,1000,2000,4000,8000,16000, 500,1000,2000,4000,8000,16000] sigma_str = "\t".join([str(int(x)) for x in Sigma])#.split(',') header_ts = 'Time [s]'+'\t'+'Normal stress [Pa]:'+'\t'+sigma_str with open(out_path + TSname + '.txt', 'w') as f: f.write(header_ts + '\n') for row in DFTS.T: csv.writer(f, delimiter='\t', lineterminator='\n').writerow(row) f.close def main(): DFTS = convert_tdms(path_in,params_conv) RST_plot(path_out, name_ts, DFTS) main()
This diff is collapsed.
 Normal stress [Pa] Shear strength [Pa] 500.0 384.19 1000.0 670.06 2000.0 1244.05 4000.0 2395.66 8000.0 4641.83 16000.0 9163.69 500.0 373.07 1000.0 667.8 2000.0 1244.02 4000.0 2387.42 8000.0 4646.12 16000.0 9112.94 500.0 377.95 1000.0 671.29