from sympy import *
import math
import numpy

a = 5.4
b = 5.4
h = 0.09
R = 225 * h
A = 1
B = 1
P_x = 0
P_y = 0
E = 2.1 * 10**5
E_1 = E
E_2 = E
my = 0.3
my_12 = my
my_21 = my
myv = (1 - my) / 2
G = E / 2 / (1 + my)
G_12 = G
G_13 = G
G_23 = G
k_x = 1 / R
k_y = 1 / R
k = 5 / 6

N = 1

for q in numpy.arange(0, 0.1, 0.1):

    x, y = symbols('x y')

    Uij = MatrixSymbol('Uij', sqrt(N), sqrt(N))
    Vij = MatrixSymbol('Vij', sqrt(N), sqrt(N))
    Wij = MatrixSymbol('Wij', sqrt(N), sqrt(N))
    Psy_xij = MatrixSymbol('Psy_xij', sqrt(N), sqrt(N))
    Psy_yij = MatrixSymbol('Psy_yij', sqrt(N), sqrt(N))

    U = 0
    V = 0
    W = 0
    Psy_x = 0
    Psy_y = 0

    for i in range(1, sqrt(N) + 1):
        for j in range(1, sqrt(N) + 1):
            X1 = sin((2 * i) * math.pi * x / a)
            X2 = sin((2 * i - 1) * math.pi * x / a)
            X3 = sin((2 * i - 1) * math.pi * x / a)
            X4 = cos((2 * i - 1) * math.pi * x / a)
            X5 = sin((2 * i - 1) * math.pi * x / a)

            Y1 = sin((2 * j - 1) * math.pi * y / b)
            Y2 = sin((2 * j) * math.pi * y / b)
            Y3 = sin((2 * j - 1) * math.pi * y / b)
            Y4 = sin((2 * j - 1) * math.pi * y / b)
            Y5 = cos((2 * j - 1) * math.pi * y / b)

            U = Uij[i-1, j-1] * X1 * Y1 + U
            V = Vij[i-1, j-1] * X2 * Y2 + V
            W = Wij[i-1, j-1] * X3 * Y3 + W
            Psy_x = Psy_xij[i-1, j-1] * X4 * Y4 + Psy_x
            Psy_y = Psy_yij[i-1, j-1] * X5 * Y5 + Psy_y

    All = BlockMatrix([[Uij], [Vij], [Wij], [Psy_xij], [Psy_yij]])
    #print(All.as_explicit())
    dUdx = diff(U, x)
    dUdy = diff(U, y)
    dVdx = diff(V, x)
    dVdy = diff(V, y)
    dWdx = diff(W, x)
    dWdy = diff(W, y)
    dPsy_xdx = diff(Psy_x, x)
    dPsy_xdy = diff(Psy_x, y)
    dPsy_ydx = diff(Psy_y, x)
    dPsy_ydy = diff(Psy_y, y)

    Teta_1 = -(1 / A * dWdx + k_x * U)
    Teta_2 = -(1 / B * dWdy + k_y * V)

    E_x = 1 / A * dUdx + 1 / A / B * V * diff(A, y) - k_x * W + 1 / 2 * Teta_1**2
    E_y = 1 / B * dVdy + 1 / A / B * V * diff(B, x) - k_y * W + 1 / 2 * Teta_2**2
    Gamma_xy = 1 / A * dVdx + 1 / B * dUdy - 1 / A / B * (U * diff(A, y) + V * diff(B, x)) + Teta_1 * Teta_2

    Hi_1 = 1 / A * dPsy_xdx + 1 / A / B * Psy_y * diff(A, y)
    Hi_2 = 1 / B * dPsy_ydy + 1 / A / B * Psy_x * diff(B, x)
    Hi_12 = 1 / 2 * (1 / A * dPsy_ydx + 1 / B * dPsy_xdy - 1 / A / B * (diff(A, y) * Psy_x + diff(B, x) * Psy_y))

    N_x = E_1 / (1 - my_12 * my_21) * (h * (E_x + my_21 * E_y))
    N_y = E_2 / (1 - my_12 * my_21) * (h * (E_y + my_12 * E_x))
    N_xy = G_12 * h * Gamma_xy
    N_yx = G_12 * h * Gamma_xy

    M_x = E_1 / (1 - my_12 * my_21) * h**3 / 12 * (Hi_1 + my_12 * Hi_2)
    M_y = E_2 / (1 - my_12 * my_21) * h**3 / 12 * (Hi_2 + my_21 * Hi_1)
    M_xy = 2 * G_12 * h**3 / 12 * Hi_12
    M_yx = 2 * G_12 * h**3 / 12 * Hi_12

    Q_x = G_13 * k * h * (Psy_x - Teta_1)
    Q_y = G_23 * k * h * (Psy_y - Teta_2)

    Esf = 1 / 2 * (N_x * E_x + N_y * E_y + 1 / 2 * (N_xy + N_yx) * Gamma_xy + M_x * Hi_1 + M_y * Hi_2 + (M_xy + M_yx) * Hi_12 + Q_x * (Psy_x - Teta_1) + Q_y * (Psy_y - Teta_2) - 2 * (q * W + P_x * U + P_y * V)) * A * B



    #Esf = 1 / 2 * (integrate(integrate((N_x * E_x + N_y * E_y + 1 / 2 * (N_xy + N_yx) * Gamma_xy + M_x * Hi_1 + M_y * Hi_2 + (M_xy + M_yx) * Hi_12 + Q_x * (Psy_x - Teta_1) + Q_y * (Psy_y - Teta_2) - 2 * (q * W + P_x * U + P_y * V)) * A * B, (x, 0, a)), (y, 0, b)))
    #print(All.as_explicit())