Commit b848a7f9 by lapm

 """ Modified on June 11, 2019 @author: Lars P. Mikkelsen, DTU Wind energy """ import numpy as np import matplotlib.pyplot as plt import warnings import sys import copy def ReadLineNumber(GeoData,SpecName,icol0=0,iskip=1,strdel=';'): # Find linenumber for SpecName in icol0, remember to use same iskip in loadtxt function later readspecs = np.loadtxt(open(GeoData,'r'),skiprows=iskip, delimiter=strdel,usecols=(icol0,),dtype='str',unpack=True) # Look for geometry data ifind=0 for j in range(len(readspecs)): if SpecName in readspecs[j]: ifind=1 iID=j break if ifind==0: sys.exit('Specimen name:'+SpecName+' does not exist in '+GeoData) return iID+iskip def FuncDelElements(x,fact=5): mask = np.ones(len(x), dtype=bool) xMax=max(x) for i in range(1,len(x)): if (100*abs(x[i]-x[i-1])/xMax > fact ) : mask[i]=False return mask def FuncCut(x,fact=5): x0=x[0] xMax=max(x) for i in range(1,len(x)): if (100*abs(x[i]-x[i-1])/xMax > fact ) : break icut=i-1 return icut def FuncRange(x,xRange): ia=0 ib=len(x) for i in range(len(x)): if x[i] > xRange[1]: ib=i-1 break for i in range(ib): if x[i] < xRange[0]: ia=i ib=ib+1 # in python the range [ia:ib] will not include ib return ia,ib def FuncBalanceStrain(x1,x2,x1Range): x2Range=(np.max(x2)*0.05,np.max(x2)*0.10) (ia2,ib2) = FuncRange(x2,x2Range) p=np.polyfit(x1[ia2:ib2] , x2[ia2:ib2] , 1) x1 = x1 + p[1]/p[0] for i in range(30): (ia1,ib1) = FuncRange(x1,x1Range) p=np.polyfit(x1[ia1:ib1] , x2[ia1:ib1], 1) x1 = x1 + p[1]/p[0] # Consider new break condition if np.abs(p[1]/p[0]) < np.max(x1)*0.0000001: break return x1 def FuncBalanceStrain2(x1,x2,x1Range): x2Range=(np.max(x2)*0.4,np.max(x2)*0.6) (ia2,ib2) = FuncRange(x2,x2Range) p=np.polyfit(x1[ia2:ib2] , x2[ia2:ib2] , 1) x1 = x1 + p[1]/p[0] for i in range(30): (ia1,ib1) = FuncRange(x1,x1Range) p=np.polyfit(x1[ia1:ib1] , x2[ia1:ib1], 1) x1 = x1 + p[1]/p[0] # Consider new break condition if np.abs(p[1]/p[0]) < np.max(x1)*0.0000001: break return x1 def TrueStress(strain,stress): stress=stress*(1.0+strain/100.0) strain=100.0*np.log(1.0+strain/100.0) return strain, stress def PlotTable(TableHeading,SpecNamesExt,TableDataExt): nrows, ncols = len(SpecNamesExt)+1, len(TableHeading) hcell, wcell = 0.5, 2. hpad, wpad = 0, 0 fig=plt.figure(figsize=((ncols+1)*wcell+wpad, nrows*hcell+hpad)) ax = fig.add_subplot(111) ax.axis('off') TableStr = [['%.4f' % j for j in i] for i in TableDataExt] #do the table the_table = ax.table(cellText=TableStr, colLabels=TableHeading,rowLabels=SpecNamesExt, loc='center') the_table.set_fontsize(24) # the_table.scale(1.5, 2) plt.plot() plt.savefig('TableOut.png', bbox_inches='tight') plt.show() def ReadSpecID(GeoData,SpecName,icol0,iskip=1,strdel=';'): # Find linenumber for SpecName in icol0, remember to use same iskip in loadtxt function later readspecs = np.loadtxt(open(GeoData,'r'),skiprows=iskip, delimiter=strdel,usecols=(icol0,),dtype='str',unpack=True) # Look for geometry data ifind=0 for j in range(len(readspecs)): if SpecName in readspecs[j]: ifind=1 iID=j break if ifind==0: sys.exit('Specimen name:'+SpecName+' does not exist in '+GeoData) return iID def ReadFloatFile(DataFile,nline,iprint=1): f=open(DataFile) lines=f.readlines() if iprint==1 : print(lines[nline]) for t in lines[nline].split(): try: a=float(t) except ValueError: pass return a def FuncStrainReplace(defl,strain,ia,ib): # to be used if error strain measure should be replace by estimate # based on fit of strain versus defl (or time) curve between index ia and ib # strain=FuncStrainReplace(defl,strain,ia,ib) p = np.polyfit(defl[ia:ib] , strain[ia:ib] ,1) for i in range(ib,len(strain)): strain[i] = defl[i]*p[0] + p[1] return strain def YieldStress02(Strain,Stress,StrainRange,YieldStrain=0.2): iyield=0 (ia,ib) = FuncRange(Strain,StrainRange) #Find stiffness pE = np.polyfit(Strain[ia:ib] , Stress[ia:ib] , 1) Emod=pE[0]/10 # Emod in GPa pEyield=np.array([0,0]) pEyield=copy(pE) #Move stiffness line with YieldStrain pEyield[1]=pE[1]-YieldStrain*pE[0] StrEmod = r'E=%4.3f GPa' % (Emod) #Find yield stress for straini,stressi in zip(Strain,Stress): if stressi < np.polyval(pEyield,straini): break iyield=iyield+1 #Save results for plot StrainY=Strain[iyield] StressY=Stress[iyield] StrStressY = r'$\sigma_Y$=%4.3f MPa' % (StressY) StrStrainY=r'$\varepsilon_Y^P$=%4.3f %%' % (YieldStrain) # Plot solution plt.figure() plt.plot(Strain,Stress,'-r',label='Data') plt.plot(Strain[ia:ib],np.polyval(pE,Strain[ia:ib]), '-g',linewidth=4,label=StrEmod) plt.plot(Strain,np.polyval(pEyield,Strain),'--k',label=StrStrainY) plt.plot(StrainY,StressY,'ko',label=StrStressY) plt.xlabel(r'$\varepsilon$ [%]',fontsize=18) plt.ylabel(r'$\sigma$ [MPa]',fontsize=18) plt.legend(loc=4) plt.xlim(np.array([-0.02,1.2])*np.max(Strain)) plt.ylim(np.array([-0.02,1.2])*np.max(Stress)) plt.grid() plt.show() return iyield def BuildTable(TableHeading,SpecNames,C1=1,C2=1,C3=1,C4=1,C5=1,C6=1,C7=1): Nspec=len(SpecNames) if C1 != 1 : C1.append(np.mean(C1[:Nspec])) C1.append(np.std(C1[:Nspec])) C1.append(np.std(C1[:Nspec])/np.sqrt(Nspec)) if C2 != 1 : C2.append(np.mean(C2[:Nspec])) C2.append(np.std(C2[:Nspec])) C2.append(np.std(C2[:Nspec])/np.sqrt(Nspec)) if C3 != 1 : C3.append(np.mean(C3[:Nspec])) C3.append(np.std(C3[:Nspec])) C3.append(np.std(C3[:Nspec])/np.sqrt(Nspec)) if C4 != 1 : C4.append(np.mean(C4[:Nspec])) C4.append(np.std(C4[:Nspec])) C4.append(np.std(C4[:Nspec])/np.sqrt(Nspec)) if C5 != 1 : C5.append(np.mean(C5[:Nspec])) C5.append(np.std(C5[:Nspec])) C5.append(np.std(C5[:Nspec])/np.sqrt(Nspec)) if C6 != 1 : C6.append(np.mean(C6[:Nspec])) C6.append(np.std(C6[:Nspec])) C6.append(np.std(C6[:Nspec])/np.sqrt(Nspec)) if C7 != 1 : C7.append(np.mean(C7[:Nspec])) C7.append(np.std(C7[:Nspec])) C7.append(np.std(C7[:Nspec])/np.sqrt(Nspec)) TableData=[] if C1 != 1 : TableData.append(C1) if C2 != 1 : TableData.append(C2) if C3 != 1 : TableData.append(C3) if C4 != 1 : TableData.append(C4) if C5 != 1 : TableData.append(C5) if C6 != 1 : TableData.append(C6) if C7 != 1 : TableData.append(C7) SpecNamesExt=[] for SpecName in SpecNames: SpecNamesExt.append(SpecName) SpecNamesExt.append('Mean') SpecNamesExt.append('STD') SpecNamesExt.append('STD of Mean') TableData=[list(i) for i in zip(*TableData)] nrows, ncols = len(SpecNamesExt)+1, len(TableHeading) hcell, wcell = 0.5, 2. hpad, wpad = 0, 0 fig=plt.figure(figsize=((ncols+1)*wcell+wpad, nrows*hcell+hpad)) ax = fig.add_subplot(111) ax.axis('off') TableStr = [['%.4f' % j for j in i] for i in TableData] #do the table the_table = ax.table(cellText=TableStr, colLabels=TableHeading,rowLabels=SpecNamesExt, loc='center') the_table.set_fontsize(18) # the_table.scale(1.1, 1.5) plt.plot() plt.savefig('TableOut.png') plt.show() \ No newline at end of file