from sympy import *
import numpy as np
import matplotlib.pyplot as plt

%matplotlib notebook
init_session()

"""
Explanation: Stresses in a disk with diametrally opposite forces
End of explanation
"""


d, P = symbols("d P", positive=True)

phi = x*atan(x/(d/2 - y)) + x*atan(x/(d/2 + y)) + P/(pi*d)*(x**2 + y**2)

Sxx = phi.diff(y, 2)
Syy = phi.diff(x, 2)
Sxy = -phi.diff(x, 1, y, 1)

display(Sxx)
display(Syy)
display(Sxy)

Sxx_fun = lambdify((x, y), Sxx.subs({P:1, d:2}), numpy)
Syy_fun = lambdify((x, y), Syy.subs({P:1, d:2}), numpy)
Sxy_fun = lambdify((x, y), Sxy.subs({P:1, d:2}), numpy)

theta, r = np.mgrid[0:2*np.pi:101j, 1e-6:1:101j]

xx = r*np.cos(theta)
yy = r*np.sin(theta)

Sxx_vec = Sxx_fun(xx, yy)
Syy_vec = Syy_fun(xx, yy)
Sxy_vec = Sxy_fun(xx, yy)

plt.figure()
plt.contourf(xx, yy, Sxx_vec, cmap="RdYlGn", vmin=-12, vmax=12)
plt.colorbar()
plt.axis("square");

plt.figure()
plt.contourf(xx, yy, Syy_vec, cmap="RdYlGn", vmin=-210, vmax=210)
plt.colorbar()
plt.axis("square");

plt.figure()
plt.contourf(xx, yy, Sxy_vec, cmap="RdYlGn")
plt.colorbar()
plt.axis("square");

"""
Explanation: The stress function for a disk of diameter $d$ with center in the origin, and radial inward and opposite forces $P$ placed at $(0, d/2)$ and $(0, -d/2)$ is given by
$$\phi = x\arctan\left[\frac{x}{d/2 - y}\right] + x\arctan\left[\frac{x}{d/2 + y}\right] + \frac{P}{\pi d}(x^2 + y^2)$$
We know that the stresses are given by
\begin{align}
\sigma_{xx} = \frac{\partial^2 \phi}{\partial x^2}\
\sigma_{yy} = \frac{\partial^2 \phi}{\partial y^2}\
\sigma_{xy} = -\frac{\partial^2 \phi}{\partial x \partial y}
\end{align}
End of explanation
"""


E, nu = symbols("E nu")
G = E/(2*(1 + nu))

Exx = simplify(1/E*(Sxx - nu*Syy))
Eyy = simplify(1/E*(Syy - nu*Sxx))
Exy = simplify(1/G*Sxy)

display(Exx)

display(Eyy)


display(Exy)

"""
Explanation: Strains
End of explanation
"""


C1, a, b = symbols("C1 a b")

ux, uy = symbols("ux uy", cls=Function)

eq1 = diff(ux(x), x) - Exx
eq1

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

"""
Explanation: Displacements
End of explanation
"""


