from IPython import display
display.SVG('KG-8.2sd.svg')

"""
Explanation: NOTE: The slope-deflection sign convention may seem strange to those used to matirx stiffness analysis, but it makes sense.  None of the slope deflection equations explicitly state a 
member 'direction' and it doesn't matter.  For example, whether you consider the column AB as
going from A to B or as going from B to A, a +ive shear at end A is still directed toward the
left. In matrix analysis, that direction matters.
Kulak & Grondin - Example 8.2
This solves a close approximation to the first-order example of Kulak and Grondin.  The major
differences are:
* the lateral resistance of the outside column stacks are ignored; only the central rigid frame is included.
* this method does not account for axial changes of length in any of the members.
The results obtained here are 
compared with those obtained via a first-order matrix analysis,
and found to agree to about 6 significant figures.
End of explanation
"""


import sympy as sy                  # use symbolic algebra
sy.init_printing()                  # print pretty math
from sdutil2 import SD, FEF         # slope-deflection library

sy.var('theta_a theta_b theta_c theta_d theta_e theta_f Delta_b Delta_c')

E = 200000
Ic = 222E6   # W310x97
Ib = 488E6   # W460x106
Hf = 21900   # Horizontal load at F
He = 43400   # Horizontal load at E

Lab = 6500   # Length of column
Lbc = 5500   # Length of column
Lbe = 10500  # Length of beam

Mab,Mba,Vab,Vba = SD(Lab,E*Ic,theta_a,theta_b,Delta_b)             # column AB, BC
Mbc,Mcb,Vbc,Vcb = SD(Lbc,E*Ic,theta_b,theta_c,Delta_c-Delta_b)

Mde,Med,Vde,Ved = SD(Lab,E*Ic,theta_d,theta_e,Delta_b)             # column DE, EF
Mef,Mfe,Vef,Vfe = SD(Lbc,E*Ic,theta_e,theta_f,Delta_c-Delta_b)

Mbe,Meb,Vbe,Veb = SD(Lbe,E*Ib,theta_b,theta_e) + FEF.udl(Lbe,55)   # beams BE, CF
Mcf,Mfc,Vcf,Vfc = SD(Lbe,E*Ib,theta_c,theta_f) + FEF.udl(Lbe,45)

eqns = [ Mba+Mbe+Mbc,      # sum of moments at B = 0
         Mcb+Mcf,          # sum of moments at C = 0
         Med+Meb+Mef,      # sum of moments at E = 0
         Mfe+Mfc,          # sum of moments at F = 0
         -Vab-Vde+Hf+He,   # sum of Fx @ base of storey 1 = 0
         -Vbc-Vef+Hf,      # sum of Fx @ base of storey 2 = 0
         theta_a,          # fixed support at A, rotation = 0
         theta_d]          # fixed support at F, rotation = 0

soln = sy.solve(eqns)
soln

"""
Explanation: Solve for joint rotations and storey translations using slope-deflection
End of explanation
"""


Mab

Mab.subs(soln)

Mab.subs(soln).n(4) * 1E-6

V = globals()   # another way to access global variables
V['Mab']

Mab is V['Mab']

V['Mab'].subs(soln).n()

"""
Explanation: Determine end moments and shears
Demonstrate how to access and convert one end moment:
End of explanation
"""


# collect the end moments in all 6 members
allm = []
V = globals()
for m in 'ab,bc,de,ef,be,cf'.split(','):
    mj = V['M'+m].subs(soln).n()*1E-6          # Mxy
    mk = V['M'+m[::-1]].subs(soln).n()*1E-6    # Myx
    allm.append((m.upper(),mj,mk))
allm

[(m,round(a,1),round(b,1)) for m,a,b in allm]  # display to one decimal place

"""
Explanation: Determine end moments by back substitution:
End of explanation
"""


# collect the end shears in all 6 members
allv = []
V = globals()
for m in 'ab,bc,de,ef,be,cf'.split(','):
    mj = V['V'+m].subs(soln).n()*1E-3          # Mxy
    mk = V['V'+m[::-1]].subs(soln).n()*1E-3    # Myx
    allv.append((m.upper(),mj,mk))
allv

[(m,round(a,1),round(b,1)) for m,a,b in allv]   # display to one decimal place

"""
Explanation: Determine end shears by back substitution:
End of explanation
"""


import pandas as pd

dd = '../matrix-methods/frame2d/data/KG82sd.d/all'   # location of mm solution

"""
Explanation: Now compare to matrix method solution (Frame2D):
End of explanation
"""


mems = pd.DataFrame(allm,columns=['ID','MZJ','MZK']).set_index('ID')
mems

"""
Explanation: Compare Moments:
Convert current member end moments solution to tabular form:
End of explanation
"""


mefs = pd.read_csv(dd+'/member_end_forces.csv').set_index('MEMBERID').loc[mems.index]
mems2 = mefs[['MZJ','MZK']] * -1E-6    # convert sign and to kN-m
mems2

"""
Explanation: Fetch solution from Frame2D:
End of explanation
"""


mdiff = (100*(1-mems/mems2))   # calculate % diff
mdiff

mdiff.abs().max()

"""
Explanation: Compare member end moments in the two solutions:
End of explanation
"""


mevs = pd.DataFrame(allv,columns=['ID','FYJ','FYK']).set_index('ID')
mevs

"""
Explanation: The maximum difference in member end moments is 0.000009% (about 7 or 8 sig figs).
Compare Shears:
Convert our end shears to tabular form:
End of explanation
"""


mevs2 = mefs[['FYJ','FYK']] * 1E-3
mevs2[['FYK']] *= -1    # change sign on end k
mevs2

"""
Explanation: Extract the end shears from Frame2D results:
End of explanation
"""


vdiff = 100*(1-mevs/mevs2)
vdiff

vdiff.abs().max()

"""
Explanation: Compare end shears in the two results:
End of explanation
"""


deltw = pd.DataFrame([('B', soln[Delta_b], soln[theta_b]),
                      ('C', soln[Delta_c], soln[theta_c])],columns=['ID','DX','RZ']).set_index('ID')
deltw

disp = pd.read_csv(dd+'/node_displacements.csv').set_index('NODEID').loc[deltw.index][['DX','RZ']]
disp['RZ'] *= -1
disp

diffd = (100*(1-deltw/disp))
diffd

diffd.abs().max()

"""
Explanation: The maximum difference is about 0.00004% which is high, but end shears on column AB are very small.
Compare Displacements
End of explanation
"""


