# coding: utf-8 from matplotlib import rcParams rcParams['font.family'] = 'sans-serif' rcParams['font.sans-serif'] = ['arial'] import matplotlib.pyplot as plt from matplotlib.patches import Polygon from mpl_toolkits.axes_grid.axislines import Subplot from math import * import sys import numpy as np def DRAWVECT(fnameF,r,z,kakom,_sigw,_angw,matno,str_t,coef): col=['#ff0000','#0000ff'] fig = plt.figure() ax1 = Subplot(fig,111) fig.add_subplot(ax1) fig.suptitle(str_t,fontsize=14,fontweight='bold') _sig=[] _rg=[] _zg=[] _ang=[] kk=0 for i in range(0,NELT): if matno[i]==1: kk=kk+1 n1=kakom[i][0]-1; p1=[r[n1],z[n1]] n2=kakom[i][1]-1; p2=[r[n2],z[n2]] n3=kakom[i][2]-1; p3=[r[n3],z[n3]] n4=kakom[i][3]-1; p4=[r[n4],z[n4]] ax1.add_patch(Polygon([p1,p2,p3,p4], closed=True,fill=True,color='#ffffcc',lw=0.05)) ax1.add_patch(Polygon([p1,p2,p3,p4], closed=True,fill=False,color='#aaaaaa',lw=0.05)) _sig=_sig+[_sigw[i]] _ang=_ang+[_angw[i]] _rg=_rg+[(r[n1]+r[n2]+r[n3]+r[n4])*0.25] _zg=_zg+[(z[n1]+z[n2]+z[n3]+z[n4])*0.25] for i in range(0,kk): if _sig[i]<0.0: icol=0 else: icol=1 al=abs(_sig[i])*coef x1=_rg[i]+0.5*al*cos(radians(_ang[i])) y1=_zg[i]-0.5*al*sin(radians(_ang[i])) x2=_rg[i]-0.5*al*cos(radians(_ang[i])) y2=_zg[i]+0.5*al*sin(radians(_ang[i])) ax1.plot([x1,x2],[y1,y2],color=col[icol],lw=0.5) txt1='Smax={0:.3f} MPa'.format(np.amax(_sig)) txt2='Smin={0:.3f} MPa'.format(np.amin(_sig)) ax1.text(0.95,0.08,txt1,fontsize=12,color='black',ha='right',va='top',transform=ax1.transAxes) ax1.text(0.95,0.04,txt2,fontsize=12,color='black',ha='right',va='top',transform=ax1.transAxes) xmin=0 xmax=20 ymin=60 ymax=100 ax1.set_xlim([xmin,xmax]) ax1.set_ylim([ymin,ymax]) ax1.set_xticks(np.arange(xmin,xmax+5,5)) ax1.set_yticks(np.arange(ymin,ymax+5,5)) ax1.set_xlabel('Distance in x-direction (m)',fontsize=12) ax1.set_ylabel('Elevation (EL.m)',fontsize=12) aspect=1.0*(ymax-ymin)/(xmax-xmin)*(ax1.get_xlim()[1] - ax1.get_xlim()[0]) / (ax1.get_ylim()[1] - ax1.get_ylim()[0]) ax1.set_aspect(aspect) ax1.axis["right"].set_visible(False) ax1.axis["top"].set_visible(False) # legend ax2 = fig.add_axes([0.38, 0.85, 0.25, 0.05],frameon=False) # X, Y, width, height ax2.plot([0,2],[1.5,1.5],color='#0000ff',lw=1) ax2.plot([0,2],[0.5,0.5],color='#ff0000',lw=1) ax2.text(3,1.5,'Tension',fontsize=10,color='black',ha='left',va='center') ax2.text(3,0.5,'Compression',fontsize=10,color='black',ha='left',va='center') ax2.set_xlim([0,10]) ax2.set_ylim([0,2]) ax2.set_xticks([]) ax2.set_yticks([]) ax2.tick_params(top="off") ax2.tick_params(bottom="off") ax2.tick_params(right="off") ax2.tick_params(left="off") # save & show the drawing plt.savefig(fnameF,dpi=300) # plt.show() plt.close() #Main routine param=sys.argv fnameR=param[1] # Data reading fin=open(fnameR,'r') text=fin.readline() text=fin.readline() nod=4 text=fin.readline() text=text.strip() text=text.split(',') NODT =int(text[0]) NELT =int(text[1]) MATEL =int(text[2]) KOZ =int(text[3]) KOR =int(text[4]) NF =int(text[5]) IPR =int(text[6]) text=fin.readline() text=fin.readline() for i in range(0,NODT): text=fin.readline() text=fin.readline() text=fin.readline() kakom = [[0]*nod for i in range(0,NELT)] for i in range(0,NELT): text=fin.readline() text=text.strip() text=text.split(',') for j in range(0,nod): kakom[i][j]=int(text[j+1]) text=fin.readline() text=fin.readline() z=[] r=[] disz=[] disr=[] for i in range(0,NODT): text=fin.readline() text=text.strip() text=text.split(',') z=z+[float(text[1])] r=r+[float(text[2])] disz=disz+[float(text[3])] disr=disr+[float(text[4])] text=fin.readline() text=fin.readline() sig1=[] sig2=[] sigt=[] angl=[] noten=[] matno=[] for i in range(0,NELT): text=fin.readline() text=text.strip() text=text.split(',') sigt=sigt+[float(text[4])] sig1=sig1+[float(text[6])] sig2=sig2+[float(text[7])] angl=angl+[float(text[8])] noten=noten+[int(text[9])] matno=matno+[int(text[10])] fin.close() # Drowing _sig1=[] _sig2=[] for i in range(0,NELT): if matno[i]==1: _sig1=_sig1+[sig1[i]] _sig2=_sig2+[sig2[i]] unit=3.0 vec1_max=max(abs(max(_sig1)),abs(min(_sig1))) vec2_max=max(abs(max(_sig2)),abs(min(_sig2))) vec_max=max(vec1_max,vec2_max) coef=unit/(vec_max*0.01) del _sig1,_sig2 fnameF='fig_vect1.png' str_t='First Principal Stress' _sigw=np.array(sig1)*0.01 _angw=np.array(angl)+90.0 DRAWVECT(fnameF,r,z,kakom,_sigw,_angw,matno,str_t,coef) fnameF='fig_vect2.png' str_t='Second Principal Stress' _sigw=np.array(sig2)*0.01 _angw=np.array(angl) DRAWVECT(fnameF,r,z,kakom,_sigw,_angw,matno,str_t,coef) sys.exit()