'Returns a Singularity Function expression which represents
the elastic curve or deflection of the Beam object.
Examples
There is a beam of length 30 meters. A moment of magnitude 120 Nm is
applied in the clockwise direction at the end of the beam. A pointload
of magnitude 8 N is applied from the top of the beam at the starting
point. There are two simple supports below the beam. One at the end
and another one at a distance of 10 meters from the start. The
deflection is restricted at both the supports.
Using the sign convention of upward forces and clockwise moment
being positive.
>>> from sympy.physics.continuum_mechanics.beam import Beam
>>> from sympy import symbols
>>> E, I = symbols(\'E, I\')
>>> R1, R2 = symbols(\'R1, R2\')
>>> b = Beam(30, E, I)
>>> b.apply_load(-8, 0, -1)
>>> b.apply_load(R1, 10, -1)
>>> b.apply_load(R2, 30, -1)
>>> b.apply_load(120, 30, -2)
>>> b.bc_deflection = [(10, 0), (30, 0)]
>>> b.solve_for_reaction_loads(R1, R2)
>>> b.deflection()
(4000*x/3 - 4*SingularityFunction(x, 0, 3)/3 + SingularityFunction(x, 10, 3)
+ 60*SingularityFunction(x, 30, 2) + SingularityFunction(x, 30, 3)/3 - 12000)/(E*I)'

def deflection(self):

    x = self.variable
    E = self.elastic_modulus
    I = self.second_moment
    if ((not self._boundary_conditions['deflection']) and (not self._boundary_conditions['slope'])):
        return ((S(1) / (E * I)) * integrate(integrate(self.bending_moment(), x), x))
    elif (not self._boundary_conditions['deflection']):
        return integrate(self.slope(), x)
    elif ((not self._boundary_conditions['slope']) and self._boundary_conditions['deflection']):
        C3 = Symbol('C3')
        C4 = Symbol('C4')
        slope_curve = (integrate(self.bending_moment(), x) + C3)
        deflection_curve = (integrate(slope_curve, x) + C4)
        bc_eqs = []
        for (position, value) in self._boundary_conditions['deflection']:
            eqs = (deflection_curve.subs(x, position) - value)
            bc_eqs.append(eqs)
        constants = list(linsolve(bc_eqs, (C3, C4)))
        deflection_curve = deflection_curve.subs({C3: constants[0][0], C4: constants[0][1]})
        return ((S(1) / (E * I)) * deflection_curve)
    C4 = Symbol('C4')
    deflection_curve = (integrate(((E * I) * self.slope()), x) + C4)
    bc_eqs = []
    for (position, value) in self._boundary_conditions['deflection']:
        eqs = (deflection_curve.subs(x, position) - value)
        bc_eqs.append(eqs)
    constants = list(linsolve(bc_eqs, C4))
    deflection_curve = deflection_curve.subs({C4: constants[0][0]})
    return ((S(1) / (E * I)) * deflection_curve)
