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


def plot_rotor(ansys):
    ansys.run("/Post1")
    ansys.run("ALLS")
    ansys.run("/erase")
    ansys.run("/win,1,ltop $/win,2,rtop")
    ansys.run("/win,3,lbot $/win,4,rbot")
    ansys.run("/win,1,off")
    ansys.run("/win,2,off")
    ansys.run("/win,3,off")
    ansys.run("/win,4,off")
    ansys.run("/win,1,on")
    ansys.plnsol("S", "X")  # Stress in Radialrichtung
    ansys.run("/noerase")
    ansys.run("/win,1,off")
    ansys.run("/win,2,on")
    ansys.run("/view,2,0,0,-1")
    ansys.plnsol("S", "X")  # Stress in Radialrichtung, Sicht von hinten
    ansys.run("/noerase")
    ansys.run("/win,2,off")
    ansys.run("/win,3,on")
    ansys.plnsol("S", "Y")  # Stress in Tangentialrichtung
    ansys.run("/noerase")
    ansys.run("/win,3,off")
    ansys.run("/win,4,on")
    ansys.plnsol("U", "SUM")
    ansys.run("/noerase")
    ansys.run("/win,4,off")
    ansys.run("/EOF")


def sim_rotor(materials, angles, outer_diameter=250, inner_diameter=50, thickness=5, omega=0.1):
    layers = len(angles)
    layer_thickness = thickness / layers

    path = os.getcwd()+"\\ws"
    mapdl = pyansys.launch_mapdl(run_location=path, interactive_plotting=True)
    # clear existing geometry
    mapdl.finish()
    mapdl.clear()

    th = 200  # [°C]        Herstellungstemperatur
    tn = 20  # [°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)
    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.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.aplot()
    mapdl.mshape(0, "2D")  # quadratic 2D Shapes
    mapdl.mshkey(1)  # 0-free meshing
    mapdl.amesh("ALL")
    # mapdl.eplot()

    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.post1()
    #mapdl.open_gui()


    # plot_rotor(mapdl)

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

    result = pyansys.read_binary(resultfile)
    # result = mapdl.result
    result.plot_principal_nodal_stress(0, 'EQV')

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

    # von Mises stress is the last column
    # must be nanmax as the shell element stress is not recorded
    maxstress = np.nanmax(stress[:, -1])

    mapdl.exit()

    return maxstress


print(sim_rotor(
    [{
        "e11": 42.5,
        "e22": 11,
        "e33": 11,
        "pr12": 0.28,
        "pr13": 0.28,
        "pr23": 0.28,
        "dens": 1950E-6,
        "g12": 4.2,
        "g13": 4.2,
        "g23": 2.56,
        "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.4E3,
        "d22": 3.1E3,
        "d33": 3.1E3
    }], [90, 0, 0, 0, 0, 0]))