from sympy import *                 # use symbolic algebra
init_printing(use_latex='mathjax')  # print nice math (use_latex arg is workaround for FormatterWarning messages)
from IPython import display

from sdutil2 import SD, FEF

"""
Explanation: Slope Deflection Example - Frame 1
Illustrates the solution of simple, statically-indeterminate framed structures
using the method of slope-deflection and symbolic algebra.
This version uses an improved version of 'sdutil' that computes end shears as well as moments.
End of explanation
"""


display.SVG('frame1.svg')

"""
Explanation: Analyze the following frame using the method of slope-deflection.
End of explanation
"""


help(SD)

# Slope deflection equations:
# theta_a,b,c,d - rotations of joints
# Delta_ab - horizontal displacement of joint b (and c) relative to a (and d)
var('EI theta_a theta_b theta_c theta_d Delta_ab')
Mab,Mba,Vab,Vba = SD(3,EI,theta_a,theta_b,Delta_ab)
Mbc,Mcb,Vbc,Vcb = SD(6,2*EI,theta_b,theta_c,0) + FEF.udl(6,36)
Mcd,Mdc,Vcd,Vdc = SD(2,EI,theta_c,theta_d,Delta_ab)

"""
Explanation: This frame has 5 independant joint displacements - rotations at the four
joints a, 
b, 
c, 
and 
d, 
and a lateral displacement of joints 
b
and 
c 
together.  Actually, we know that joint d does not rotate, but we
will pretend it does for now, and will handle that as a displacement constraint
later.
So now we write the slope deflection equations.  These express the end moments
as functions of the five unknown displacements.  Note also that we do not treat the
pin-ended member ab specially -- it is not necessary to do so.
We use a very simple utility function SD() to compute the end moments and shears of a segment,
as a function of the end displacements:
End of explanation
"""


Mbc

"""
Explanation: Display just one of the end moments to see what the values look like.  They are indeed algebraic expressions.
End of explanation
"""


eqns = [Mab, Mba+Mbc, Mcb+Mcd, Vab+Vdc-60, theta_d]
soln = solve(eqns,[theta_a,theta_b,theta_c,theta_d,Delta_ab])
soln

"""
Explanation: Now we write five linear equations to solve for the five unknown displacements.
Three are equilibrium equations and involve moments at three of the joints: $Mab = 0$, $Mba+Mbc=0$, and $Mcb+Mcd=0$.
One equation involves the equilibrium of the column shears with the horizontal applied
load: $Vab+Vdc-60=0$.  The fifth equation is just the known rotation at joint d: $\theta_d=0$.
In sympy, we build a list of the five equations, all expressed as $f(..,..) = 0$ and giving only the left side of each (the right side always being 0 is not given).  We then call the solver and give it a list of variables to solve for:
End of explanation
"""


{v:eval(v).subs(soln).n(4) for v in 'Mab,Mba,Mbc,Mcb,Mcd,Mdc,Vab,Vdc'.split(',')}

"""
Explanation: The above solution is returned as a Python dictionary, making it easy to substitute for those
variables in other equations.
In fact, that is done now, substituting the displacement solutions into each
of the slope-deflection equations.  The result is a list of the member end moments (and shears).
The '.n(4)' method forces a numerical answer with 4 significant figures.
End of explanation
"""


Ra = Vbc
Rd = -Vcb
Ha = -Vab
Hd = -Vdc
Md = Mdc
[v.subs(soln).n(4) for v in [Ra,Ha,Rd,Hd,Md]]

"""
Explanation: And now - the reactions at supports a and d:
End of explanation
"""


# sum horizontal forces
(Ha + Hd + 60).subs(soln)

# sum vertical forces
(Ra + Rd - 36*6).subs(soln)

# sum moments about a
(60*3 + 36*6*3 + Hd*1 + Md - Rd*6).subs(soln)

"""
Explanation: Equilibrium checks:
End of explanation
"""


%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

var('x')
Mbcx = (Ra*x + Vab*3 - 36*x*x/2).subs(soln)  # M as a function of x
fn = lambdify((x),Mbcx)   # convert moment expression to a function of one argument
z = np.linspace(0,6)
plt.plot(z,fn(z),[0,6],[0,0])
plt.fill_between(z,fn(z),color='lightgrey')
##plt.show()

"""
Explanation: Plots
Plot bending moment in bc as a function of distance x from left end, b:
End of explanation
"""


[r.n(4) for r in solve(Mbcx)]

"""
Explanation: Find the locations along span bc of zero moment:
End of explanation
"""


xs = solve([Mbcx.diff(x)])
xs

xs[x].n(4), Mbcx.subs(xs).n(4)

"""
Explanation: Find point of maximum moment and the maximum moment in the span:
End of explanation
"""


