from scipy.interpolate import interp1d
import matplotlib.path as mpltPath
from FEM.PlaneStress import PlaneStress
from FEM.Mesh.Delaunay import Delaunay
import numpy as np
import matplotlib.pyplot as plt
from funciones import empujes
import copy
h = 0.5
Kh = 1E5

L = 15
H = 0.6
# coords = [[0, 0], [L, 0], [L, H], [0, H]]

params = Delaunay._strdelaunay(True, True, '0.01', o=1)

# cb = geometria.cbFromSegment()
# cbn = [[0, 0.0], [1, 0.0], [-1, 0.0], [-2, 0.0]]
extra_x = 7.0
extra_y = 5.0
w_d = H+2.0*extra_x
h_d = L+extra_y
verticesTotal = [[0.0, 0.0], [w_d, 0.0], [w_d, h_d], [0.0, h_d]]
vertextra = [[extra_x, h_d], [extra_x, h_d-L],
             [extra_x+H, h_d-L], [extra_x+H, h_d]]
segextra = [[0, 1], [1, 2], [2, 3], [3, 0]]

holes = []
# If holes center is outside the domains, the hole region will be highlighted
hole = {'center': [-100, -100], 'segments': segextra, 'vertices': vertextra}
holes += [hole]

params = Delaunay._strdelaunay(constrained=True, delaunay=True, a='0.1', o=2)

geometria = Delaunay(verticesTotal, params, nvn=2, holes_dict=holes)
cb = geometria.cbFromSegment(3, 0.0, nv=1)
cb += geometria.cbFromSegment(1, 0.0, nv=1)
cb += geometria.cbFromSegment(0, 0.0, nv=2)
geometria.setCbe(cb)
E = 27800000
v = 0.21
E_suelo = Kh
v_suelo = 0.5
EE = []
VV = []
for centroide in geometria.centroids:
    path = mpltPath.Path(vertextra)
    inside2 = path.contains_points([centroide])
    EE += [E]
    VV += [v]
    if inside2[0]:
        EE[-1] = E_suelo
        VV[-1] = v_suelo

# geometria.show()
# plt.show()


err = 1.0
def UWX(xx): return 0.0


U1p1 = 0.0
for i in range(10):
    z, _, _, _, _, pw, P0, Pa, Pp = empujes(
        15, n=500, h=h, NF=[3, 3], q_sobrecarga=0.0)
    pwI = interp1d(z, pw[1:, 0])
    pwD = interp1d(z, pw[1:, 1])
    P0I = interp1d(z, P0[1:, 0])
    P0D = interp1d(z, P0[1:, 1])
    PAI = interp1d(z, Pa[1:, 0])
    PAD = interp1d(z, Pa[1:, 1])
    PPI = interp1d(z, Pp[1:, 0])
    PPD = interp1d(z, Pp[1:, 1])

    def W(x):
        PI = 0
        PD = 0
        ux = UWX(x)
        if ux >= 0:
            sigmaII = P0I(x)-Kh*ux
        else:
            sigmaII = P0I(x)+Kh*ux
        if ux < 0:
            sigmaDD = (P0D(x)+Kh*ux)*(x > h)
        else:
            sigmaDD = (P0D(x)-Kh*ux)*(x > h)

        if ux >= 0:
            if sigmaII < PAI(x):
                sigmaII = PAI(x)
                PI = 1
            if sigmaDD > PPD(x):
                sigmaDD = PPD(x)
                PD = 1
        else:
            if sigmaII > PPI(x):
                sigmaII = PPI(x)
                PI = 1
            if sigmaDD < PAD(x):
                sigmaDD = PAD(x)
                PD = 1
        return pwI(x) - pwD(x) + P0I(x)*(1-PI) + sigmaII*PI - P0D(x)*(1-PD) - sigmaDD*(PD)
    geometria.loadOnSegment(4, fy=lambda x: W(x))
    geometria.show()
    plt.show()

    def cf(x):
        PI = 0
        PD = 0
        ux = UWX(x)
        if ux >= 0:
            sigmaII = P0I(x)-Kh*ux
        else:
            sigmaII = P0I(x)+Kh*ux
        if ux < 0:
            sigmaDD = (P0D(x)+Kh*ux)*(x > h)
        else:
            sigmaDD = (P0D(x)-Kh*ux)*(x > h)

        if ux >= 0:
            if sigmaII < PAI(x):
                PI = 1
            if sigmaDD > PPD(x):
                PD = 1
        else:
            if sigmaII > PPI(x):
                PI = 1
            if sigmaDD < PAD(x):
                PD = 1
        if PD or PI:
            print(x, PD, PI, Kh*(1-PI)+Kh*(x >= h)*(1-PD), pwI(x) - pwD(x) +
                  P0I(x)*(1-PI) + sigmaII*PI - P0D(x)*(1-PD) - sigmaDD*(PD))
        return Kh*(1-PI)+Kh*(x >= h)*(1-PD)
    # original_nodes = geometria.giveNodesOfSegment(0, 10**(-10))
    # coords_gdl_convective = []
    # for aaaaa in original_nodes:
    #     coords_gdl_convective += [geometria.gdls[aaaaa][0]]
    # gdl_convective = original_nodes*geometria.nvn+(2-1)

    O = PlaneStress(geometria, EE, VV, 1.0, verbose=True)
    O.elementMatrices()
    O.ensembling()
    O.borderConditions()
    O.solveES()
    # X, U1, U2, U3, U4 = O.postProcess(plot=True)
    X, U1 = O.profile([extra_x+H/2, h_d-L], [extra_x+H/2, h_d], n=150)
    plt.show()
    U1 = np.array(U1)[:, 0, 0]
    UWX = interp1d(X, U1)
    err = np.max(np.abs(np.array(U1)-np.array(U1p1)))
    U1p1 = copy.deepcopy(U1)
    if err < 1*10**(-10):
        plt.show()
        break

"""
THIS SHIT DOES NOT WORKS."""
