#***************************************
# make stress vector diagram for FEM
#***************************************
# 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 DRAWVEC(fnameF,x,y,kakom,_sigw,_ang,str_t,coef):
    col=['#ff0000','#0000ff']
    fig = plt.figure()
    ax1 = Subplot(fig,111)
    fig.add_subplot(ax1)
    xg=np.zeros(NELT)
    yg=np.zeros(NELT)
    for i in range(0,NELT):
        n1=kakom[i][0]-1
        n2=kakom[i][1]-1
        n3=kakom[i][2]-1
        n4=kakom[i][3]-1
        p1=[x[n1],y[n1]]
        p2=[x[n2],y[n2]]
        p3=[x[n3],y[n3]]
        p4=[x[n4],y[n4]]
        ax1.add_patch(Polygon([p1,p2,p3,p4], closed=True,fill=True,color='#ffffcc',lw=0.1))
        ax1.add_patch(Polygon([p1,p2,p3,p4], closed=True,fill=False,color='#aaaaaa',lw=0.1))
        xg[i]=np.average([x[n1],x[n2],x[n3],x[n4]])
        yg[i]=np.average([y[n1],y[n2],y[n3],y[n4]])
    for i in range(0,NELT):
        if _sigw[i]<0.0:
            icol=0
        else:
            icol=1
        al=abs(_sigw[i])*coef
        x1=xg[i]-0.5*al*cos(radians(_ang[i]))
        y1=yg[i]-0.5*al*sin(radians(_ang[i]))
        x2=xg[i]+0.5*al*cos(radians(_ang[i]))
        y2=yg[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(_sigw))
    txt2='Smin={0:.3f} MPa'.format(np.amin(_sigw))
    ax1.text(0.95,0.95,txt1,fontsize=12,color='black',ha='right',va='top',transform=ax1.transAxes)
    ax1.text(0.95,0.90,txt2,fontsize=12,color='black',ha='right',va='top',transform=ax1.transAxes)
    xmin=-12
    xmax=12
    ymin=62
    ymax=82
    ax1.set_xlim([xmin,xmax])
    ax1.set_ylim([ymin,ymax])
    ax1.set_xticks(np.arange(xmin,xmax+2,2))
    ax1.set_yticks(np.arange(ymin,ymax+2,2))
    ax1.set_xlabel('Distance in x-direction (m)',fontsize=12)
    ax1.set_ylabel('Elevation (EL.m)',fontsize=12)
    aspect = (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)
    ax1.set_title(str_t)
    # legend
    ax2 = fig.add_axes([0.2, 0.8, 0.2, 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=200)
#    plt.show()
    plt.close()


#Main routine
param=sys.argv
fnameR=param[1]
#fnameW=param[2]

# Data input
fin=open(fnameR,'r')
text=fin.readline()
text=fin.readline()

text=fin.readline()
text=text.strip()
text=text.split(',')
nod    =int(text[0])
NODT   =int(text[1])
NELT   =int(text[2])
MATEL  =int(text[3])
KOX    =int(text[4])
KOY    =int(text[5])
NF     =int(text[6])
NSTRESS=int(text[7])
IPR    =int(text[8])

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()
x=[]
y=[]
disx=[]
disy=[]
for i in range(0,NODT):
    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    x=x+[float(text[1])]
    y=y+[float(text[2])]
    disx=disx+[float(text[3])]
    disy=disy+[float(text[4])]
text=fin.readline()
text=fin.readline()
sig1=[]
sig2=[]
angl=[]
noten=[]
matno=[]
for i in range(0,NELT):
    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    sig1=sig1+[float(text[7])]
    sig2=sig2+[float(text[8])]
    angl=angl+[float(text[9])]
    noten=noten+[int(text[10])]
    matno=matno+[int(text[11])]
fin.close()

# Drowing
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)

fnameF='fig_vect1.png'
str_t='First Principal Stress Vector'
_sigw=np.array(sig1)*0.01
_ang=np.array(angl)
DRAWVEC(fnameF,x,y,kakom,_sigw,_ang,str_t,coef)

fnameF='fig_vect2.png'
str_t='Second Principal Stress Vector'
_sigw=np.array(sig2)*0.01
_ang=np.array(angl)+90.0
DRAWVEC(fnameF,x,y,kakom,_sigw,_ang,str_t,coef)