# https://akaszynski.github.io/pyansys/ansys_control.html#initial-setup-and-example
# pip install pyansys --upgrade
import pyansys
import numpy as np
import os
import time


def sim_rotor(materials, angles, outer_diameter=0.250, inner_diameter=0.050, thickness=0.005, omega=1):
    layers = len(angles)
    layer_thickness = thickness / layers

    path = os.getcwd() + "\\ws"
    mapdl = pyansys.launch_mapdl(run_location=path, interactive_plotting=True, loglevel='ERROR', override=True)

    # clear existing geometry
    mapdl.finish()
    mapdl.clear()

    th = 50  # [°C]        Herstellungstemperatur
    tn = 50  # [°C]        Nutzungstemperatur

    mapdl.prep7()
    # define Material
    for i in range(0, len(materials)):
        mapdl.mptemp(i + 1, th)

        mapdl.mpdata("EX", i + 1, "", materials[i]["e11"])
        mapdl.mpdata("EY", i + 1, "", materials[i]["e22"])
        mapdl.mpdata("EZ", i + 1, "", materials[i]["e33"])

        mapdl.mpdata("PRXY", i + 1, "", materials[i]["pr12"])
        mapdl.mpdata("PRXZ", i + 1, "", materials[i]["pr13"])
        mapdl.mpdata("PRYZ", i + 1, "", materials[i]["pr23"])

        mapdl.mpdata("DENS", i + 1, "", materials[i]["dens"])

        mapdl.mpdata("GXY", i + 1, "", materials[i]["g12"])
        mapdl.mpdata("GXZ", i + 1, "", materials[i]["g13"])
        mapdl.mpdata("GYZ", i + 1, "", materials[i]["g23"])

        # Waermeausdenungsk
        mapdl.mpdata("ALPX", i + 1, "", materials[i]["alp11"])
        mapdl.mpdata("ALPY", i + 1, "", materials[i]["alp22"])
        mapdl.mpdata("ALPZ", i + 1, "", materials[i]["alp33"])

        # Waermeleitfaehigk
        mapdl.mpdata("KXX", i + 1, "", materials[i]["k11"])
        mapdl.mpdata("KYY", i + 1, "", materials[i]["k22"])
        mapdl.mpdata("KZZ", i + 1, "", materials[i]["k33"])

        # Quellkoeffizient
        mapdl.mpdata("BETX", i + 1, "", materials[i]["bet11"])
        mapdl.mpdata("BETY", i + 1, "", materials[i]["bet22"])
        mapdl.mpdata("BETZ", i + 1, "", materials[i]["bet33"])

        # Diffusionskoeffizient
        mapdl.mpdata("DXX", i + 1, "", materials[i]["d11"])
        mapdl.mpdata("DYY", i + 1, "", materials[i]["d22"])
        mapdl.mpdata("DZZ", i + 1, "", materials[i]["d33"])

    # define Element
    mapdl.et(1, "SHELL181")
    # enable Laminate
    mapdl.run("KEYOPT,1,8,1")

    # define Laminate
    mapdl.run("SECT, 1, SHELL")

    # secdata, Thickness ,MatId,theta,Number of integration points in layer
    if len(angles) != len(materials):
        for i in range(0, len(angles)):
            mapdl.secdata(layer_thickness, 1, angles[i], 3)
            print(str(angles[i]))
    else:
        for i in range(0, len(angles)):
            mapdl.secdata(layer_thickness, i + 1, angles[i], 3)

    # Shell node will be offset to midplane of the section
    mapdl.secoffset("MID")
    mapdl.run("seccontrol,,,,,,,")

    mapdl.run("/ESHAPE,1.0")
    # define Geometry
    # change global coordinate system to cylindric system
    mapdl.csys(1)
    mapdl.cyl4(0, 0, inner_diameter / 2, 0, outer_diameter / 2, 180)
    mapdl.cyl4(0, 0, inner_diameter / 2, 180, outer_diameter / 2, 360)
    # mapdl.aplot()
    mapdl.nsel("S", "LOC", "Y", 0)
    mapdl.nummrg("ALL")
    mapdl.nsel("S", "LOC", "Y", 180)
    mapdl.nummrg("ALL")
    mapdl.lsel("S", "", "", 1, 7, 2)
    mapdl.lesize("ALL", "", "", 20)
    mapdl.lsel("S", "", "", 2, 4, 2)
    mapdl.lesize("ALL", "", "", 20)
    mapdl.lsel("ALL")
    mapdl.mshape(0, "2D")  # quadratic 2D Shapes
    mapdl.mshkey(1)  # 0-free meshing
    mapdl.amesh("ALL")

    mapdl.run("/SOLU")
    mapdl.omega("", "", omega)  # [rad/s]	Winkelgeschwindigkeit(+: gegen Uhrzeigersinn))
    mapdl.nsel("S", "LOC", "X", inner_diameter / 2)
    mapdl.d("ALL", "UX", 0)
    mapdl.d("ALL", "UY", 0)
    mapdl.d("ALL", "UZ", 0)
    mapdl.d("ALL", "ROTX", 0)
    mapdl.d("ALL", "ROTY", 0)
    mapdl.d("ALL", "ROTZ", 0)
    mapdl.run("ALLS")

    mapdl.outres("ALL", "ALL")
    # mapdl.open_gui()
    mapdl.solve()
    mapdl.finish()

    mapdl.run("/Post1")

    result = mapdl.result
    mapdl.aplot()
    # result.plot_principal_nodal_stress(0, '1', show_edges=True, show_axes=True)
    result.plot_nodal_solution(0, 'x')
    result.plot_principal_nodal_stress(0, '1', show_edges=True, show_axes=True)
    # """
    nodenump, stress = result.principal_nodal_stress(0)

    # von first principle stress is the first column
    # must be nanmax as the shell element stress is not recorded


    mapdl.run('*GET, NODE_COUNT, NODE, ,COUNT')
    mapdl.run('*GET, ELEM_COUNT, ELEM, ,COUNT')
    mapdl.load_parameters()
    print(mapdl.parameters['NODE_COUNT'])
    print(mapdl.parameters['ELEM_COUNT'])

    s1max = angles
    s2max = angles
    srmax = angles
    stmax = angles

    for i in range(0, len(angles)):
        mapdl.run("LAYER, " + str(i + 1))  #
        mapdl.run("*DIM, DIM_TEST, ARRAY, ELEM_COUNT, 2")
        mapdl.run("ETABLE, TEST, S, 1, MAX")
        mapdl.run("*VGET, DIM_TEST(1,1), ETAB, TEST")

        mapdl.run("*DIM, OUT_STRESS, ARRAY, NODE_COUNT, 4")
        mapdl.run("*VGET, OUT_STRESS(1, 1), NODE,, S, 1")
        mapdl.run("*VGET, OUT_STRESS(1, 2), NODE,, S, 2")
        mapdl.run("*VGET, OUT_STRESS(1, 3), NODE,, S, X")
        mapdl.run("*VGET, OUT_STRESS(1, 4), NODE,, S, Y")
        params, arrays = mapdl.load_parameters()
        print(arrays['TEST'][0])
        s1max[i] = np.nanmax(arrays['OUT_STRESS'][:, 0])
        s2max[i] = np.nanmax(arrays['OUT_STRESS'][:, 1])
        srmax[i] = np.nanmax(arrays['OUT_STRESS'][:, 2])
        stmax[i] = np.nanmax(arrays['OUT_STRESS'][:, 3])

        # print(str(element_stress[0]))
        # s1max = np.nanmax(stresses)
        # nodenump, stress = result.principal_nodal_stress(0)
        # s1max[i] = np.nanmax(element_stress[:, -5])

    s1max = np.nanmax(s1max)
    s2max = np.nanmax(s2max)
    srmax = np.nanmax(srmax)
    stmax = np.nanmax(stmax)

    # von second principle stress is the first column
    # must be nanmax as the shell element stress is not recorded
    #s2max = np.nanmax(stress[:, -4])

    #nodenum, stress = result.nodal_stress(0)
    # von Mises stress is the last column
    # must be nanmax as the shell element stress is not recorded
    #srmax = np.nanmax(stress[:, -6])
    #stmax = np.nanmax(stress[:, -5])

    print(str(s1max))

    # """

    """
    # access results using ANSYS object
    resultfile = os.getcwd() + '\\ws\\file.rst'
    result = pyansys.read_binary(resultfile)

    # plot_nodal_solution doesnt work
    nodenump, stress = result.principal_nodal_stress(0)

    # von first principle stress is the first column
    # must be nanmax as the shell element stress is not recorded
    s1max = np.nanmax(stress[:, -5])

    # von second principle stress is the first column
    # must be nanmax as the shell element stress is not recorded
    s2max = np.nanmax(stress[:, -4])

    nodenum, stress = result.nodal_stress(0)
    # von Mises stress is the last column
    # must be nanmax as the shell element stress is not recorded
    srmax = np.nanmax(stress[:, -6])
    stmax = np.nanmax(stress[:, -5])

    print(str(s1max))
    # result.plot_nodal_solution(0)
    # """
    mapdl.exit()
    time.sleep(0.5)

    return s1max, s2max, srmax, stmax


print(sim_rotor(
    [{
        "e11":      42.5E9  ,
        "e22":      11E9    ,
        "e33":      11E9    ,
        "pr12":     0.28    ,
        "pr13":     0.28    ,
        "pr23":     0.28    ,
        "dens":     1950    ,
        "g12":      4.2E9   ,
        "g13":      4.2E9   ,
        "g23":      2.56E9  ,
        "alp11":    5.7E-6  ,
        "alp22":    45E-6   ,
        "alp33":    45E-6   ,
        "k11":      0.72E3  ,
        "k22":      0.5E3   ,
        "k33":      0.5E3   ,
        "bet11":    0E-3,
        "bet22":    4E-3,
        "bet33":    4E-3,
        "d11":      4.4E-3,
        "d22":      3.1E-3,
        "d33":      3.1E-3
    }], [90, 0, 0, 0, 0, 0]))