from __future__ import division
import math
import numpy as np
from sympy import solve, symbols, pprint
from sympy import diff
x,y=symbols('x,y')
E=200000 #youngs modulus (N/mm^2)
nu=0.3 #poissons ratio
u1=0.001
u2=0.003
u3=-0.003 #displacements of corners (m)
u4=0
v1=-0.004
v2=-0.002
v3=0.001
v4=0.001
a=np.array([[1,-2,-1,2],[1,2,-1,-2],[1,2,1,2],[1,-2,1,-2]])
b=np.array([u1,u2,u3,u4])
alpha=np.linalg.solve(a,b)
alpha1=alpha[0]
alpha2=alpha[1]
alpha3=alpha[2]
alpha4=alpha[3]
a=np.array([[1,-2,-1,2],[1,2,-1,-2],[1,2,1,2],[1,-2,1,-2]])
b=np.array([v1,v2,v3,v4])
alpha=np.linalg.solve(a,b)
alpha5=alpha[0]
alpha6=alpha[1]
alpha7=alpha[2]
alpha8=alpha[3]
u=alpha1+alpha2*x+alpha3*y+alpha4*x*y
v=alpha5+alpha6*x+alpha7*y+alpha8*x*y
ex=diff(u,x)
ey=diff(v,y)
Yxy=diff(u,y)+diff(v,x)
ex=-0.000125
ey=.002
Yxy=-0.0015
sigma_x=E/(1-nu**2)*(ex+nu*ey)
sigma_y=E/(1-nu**2)*(ey+nu*ex)
print "longitudinal stress in x direction = %3.1f N/mm^2"%(sigma_x)
print "longitudinal stress in y direction = %3.1f N/mm^2"%(sigma_y)
T_xy=E/2/(1+nu)*Yxy
print "shear stress at the center of plate = %3.1f N/mm^2"%(T_xy)