import math
import sys
assert mw.version() == 13
nodes = [] # create an empty list for all nodes
lx = [] # create an empty list for all x nodes
ly = [] # create an empty list for all y nodes
lz = [] # create an empty list for all z nodes
m_u = [] # create an empty list for all Mu bending moment
#mw.message(str(len(mw.selected_nodes())))
if len(mw.selected_nodes()) < 3: # check if nodes have been selected
    sys.exit("You need at least 3 nodes to be selected!")
g = mw.input("Input Bending Moment type (only letter 'u' or 'v'):") # input the type of moment: on U or V direction
if g == "u":
    g = "momentu"
else:
    g = "momentv"
for j in mw.selected_nodes(): # get selected nodes from Mecway     
    x = mw.node(j).x # get x value from every node
    y = mw.node(j).y # get y value from every node
    z = mw.node(j).z # get z value from every node
    m_u.append(solution.node_value(str(g), j)) # gets "momentu" (bending moment/length on U local axis) in each node in 'Newton * meter'
    lx.append(x) # add the x values to the empty list
    ly.append(y) # add the y values to the empty list
    lz.append(z) # add the z values to the empty list
dist = math.sqrt((max(lx) - min(lx))**2 + (max(ly) - min(ly))**2 + (max(lz) - min(lz))**2) # distance between 2 points in cartesian space
#m_u = [num for num in m_u if num >= 0] # remove negative values of mu because they represent compression
mu = abs(sum(m_u) / len(m_u)) # get the avarage value for Mu. Is absolute, in case of negative moments.  
for el in solution.all_elements():
    for ii in solution.nodes(el):
        if j == ii:
            id = int(el) # convert string into integer
            density = mw.mass(id)/mw.volume(id) # calculate density of the element
            surf = mw.surface_area(FaceId(id, 5)) # get surface area 
            thk = mw.volume(id)/mw.surface_area(FaceId(id, 5))
# Calculate Bending Moment Capacity in 'Netown' and 'meter':
#Input data
reo_cover = 50 #mm - cover to reinforcement
bmr_diam = 12 #mm - diameter of the reo bar
bmr_phi = 0.8 # Bending Moment capacity reduction coefficient
bmr_Fsy = 500 #MPa -  reinforcement steel strength in MPa
bmr_Fc = 50 #compression concrete strength in MPa
#Output data
reo_cover = float(reo_cover) / 1000 # convert in meters
bmr_diam = float(bmr_diam) / 1000 # convert in meters
bmr_Fsy = float(bmr_Fsy) * 1000000
bmr_Fc = float(bmr_Fc) * 1000000
bmr_long_sect = dist #m theoretical width of the section
bmr_eff_depth = thk - reo_cover - float(bmr_diam / 2) # effective depth
BMR_req = (1 / float(bmr_phi)) * mu # Adjust Mu with the invert of reduction factor coefficient to get Bending Moment Capacity
# Calculate area of reinforcement required for Mu
x_req = bmr_eff_depth * (1 - math.sqrt(1 - (2 * BMR_req) / (bmr_long_sect * (bmr_eff_depth**2) * bmr_Fc)))
bmr_A_st = BMR_req / (bmr_Fsy * (bmr_eff_depth - (x_req / 2))) # area of reinforcement corresponding to the moment
bmr_no_bars = int(math.ceil(4 * (bmr_A_st * 1000000) / (math.pi * (bmr_diam * 1000)**2))) # get the number of bars
# Re-calculate the Bending Moment Capacity of the section, based on the exact number of bars
bmr_A_st = bmr_no_bars * (math.pi * ((bmr_diam * 1000)**2) / 4) * 0.0000010 # convert from mm^2 in meter^2
BMR_req2 = bmr_phi * bmr_A_st * bmr_Fsy * bmr_eff_depth * (1 - (0.6 * bmr_A_st * bmr_Fsy)/(bmr_long_sect * bmr_eff_depth * bmr_Fc)) # Bending Moment Capacity of the section in N * m
#Show final message
mw.message("Min. requirement: " + str(bmr_no_bars) + " x N" + str(int(bmr_diam * 1000)) + " bars for an avarage BM of: " + str(round(float(mu/1000),2)) + " kNm for section: " + str(int(dist*1000)) + " x " + str(int(1000*thk)) + " mm^2, for a BMC of: " + str(round(float(BMR_req2/1000),2)) + "kNm.")