from sympy import *

%matplotlib notebook
init_printing()

x = symbols('x')
E, I = symbols('E I', positive=True)
C1, C2, C3, C4 = symbols('C1 C2 C3 C4')
w, M, q, f = symbols('w M q f', cls=Function)
EI = symbols('EI', cls=Function, nonnegative=True)


M_eq = -diff(M(x), x, 2) - q(x)

M_eq

M_sol = dsolve(M_eq, M(x)).rhs.subs([(C1, C3), (C2, C4)])

M_sol

w_eq = f(x) + diff(w(x),x,2)
w_eq

w_sol = dsolve(w_eq, w(x)).subs(f(x), M_sol/EI(x)).rhs

w_sol

"""
Explanation: Euler Bernoulli Beam "solver"
The Euler-Bernoulli equation describes the relationship between the beam's deflection and the applied load
$$\frac{d^2}{dx^2}\left(EI\frac{d^2w}{dx^2}\right) = q \enspace .$$
The curve $w(x)$ describes the delection of the beam at some point $x$, $q$ is a distributed load.
This equation cannot be solve in this form in Sympy. Nevertheless, we can "trick" it to do it for us. Let us rewrite the equation as two equations
$$\begin{align}
&-\frac{d^2 M}{dx^2} = q \enspace ,\
&- \frac{d^2w}{dx^2} = \frac{M}{EI} \enspace ,
\end{align}$$
where $M$ is the bending moment in the beam. We can, then, solve the two equation as if they have source terms and then couple the two solutions.
The code below do this
End of explanation
"""


sub_list = [(q(x), 0), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()

L, F = symbols('L F')
# Fixed end
bc_eq1 = w_sol1.subs(x, 0)
bc_eq2 = diff(w_sol1, x).subs(x, 0)
# Free end
bc_eq3 = diff(w_sol1, x, 2).subs(x, L)
bc_eq4 = diff(w_sol1, x, 3).subs(x, L) + F/(E*I)

[bc_eq1, bc_eq2, bc_eq3, bc_eq4]

constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])
constants

w_sol1.subs(constants).simplify()

"""
Explanation: We want to be sure that this solution is ok. We replaced known values for $E$, $I$ and $q$ to check it.
Cantilever beam with end load
End of explanation
"""


sub_list = [(q(x), 1), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()

L = symbols('L')
# Fixed end
bc_eq1 = w_sol1.subs(x, 0)
bc_eq2 = diff(w_sol1, x).subs(x, 0)
# Free end
bc_eq3 = diff(w_sol1, x, 2).subs(x, L)
bc_eq4 = diff(w_sol1, x, 3).subs(x, L)

constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])

w_sol1.subs(constants).simplify()

"""
Explanation: Cantilever beam with uniformly distributed load
End of explanation
"""


sub_list = [(q(x), exp(x)), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()

L = symbols('L')
# Fixed end
bc_eq1 = w_sol1.subs(x, 0)
bc_eq2 = diff(w_sol1, x).subs(x, 0)
# Free end
bc_eq3 = diff(w_sol1, x, 2).subs(x, L)
bc_eq4 = diff(w_sol1, x, 3).subs(x, L)

constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])

w_sol1.subs(constants).simplify()

"""
Explanation: Cantilever beam with exponential loading
End of explanation
"""


k = symbols('k', integer=True)
C = symbols('C1:4')
D = symbols('D', cls=Function)

w_sol1 = 6*(C1 + C2*x) - 1/(E*I)*(3*C3*x**2 + C4*x**3 -
                             6*Sum(D(k)*x**(k + 4)/((k + 1)*(k + 2)*(k + 3)*(k + 4)),(k, 0, oo)))

w_sol1

"""
Explanation: Load written as a Taylor series and constant EI
We can prove that the general function is written as
End of explanation
"""


Q, alpha = symbols("Q alpha")
sub_list = [(q(x), Q), (EI(x), E*x**3/12/tan(alpha))]
w_sol1 = w_sol.subs(sub_list).doit()

M_eq = -diff(M(x), x, 2) - Q

M_eq

M_sol = dsolve(M_eq, M(x)).rhs.subs([(C1, C3), (C2, C4)])

M_sol

w_eq = f(x) + diff(w(x),x,2)
w_eq

w_sol1 = dsolve(w_eq, w(x)).subs(f(x), M_sol/(E*x**3/tan(alpha)**3)).rhs

w_sol1 = w_sol1.doit()

expand(w_sol1)

limit(w_sol1, x, 0)

L = symbols('L')
# Fixed end
bc_eq1 = w_sol1.subs(x, L)
bc_eq2 = diff(w_sol1, x).subs(x, L)
# Finite solution
bc_eq3 = C3

constants = solve([bc_eq1, bc_eq2, bc_eq3], [C1, C2, C3, C4])

simplify(w_sol1.subs(constants).subs(C4, 0))

"""
Explanation: Uniform load and varying cross-section
End of explanation
"""


M = -E*x**3/tan(alpha)**3*diff(w_sol1.subs(constants).subs(C4, 0), x, 2)
M

diff(M, x)

w_plot = w_sol1.subs(constants).subs({C4: 0, L: 1, Q: -1, E: 1, alpha: pi/9})
plot(w_plot, (x, 1e-6, 1));

from IPython.core.display import HTML
def css_styling():
    styles = open('./styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()

"""
Explanation: The shear stress would be
End of explanation
"""


