Week 11: unset_show bug and documentation

Last week, I created #14967 for implementation of plotting methods. Soon after pushing my commits, many of the jobs failed on Travis. It was strange as I was not able to reciprocate the errors locally.

After discussing it on Gitter, I got to know that it was due to the printing of plots using TextBackend in the doctest in absence of matplotlib. As matplotlib was present in my system,  doctest used matplotlib backend instead of TextBackend locally, hence passing all tests. Kalevi suggested using unset_show to stop the printing of plots during doctest but apparently, unset_show didn’t work for TextBackend. This was fixed by #14984 later that day and #14967 passed all the tests after former one was merged.

This week, I also started editing #14453 for documentation. It included a few beam problems along with their ascii diagrams.

Next Week

  • Make sure #14967 and#14453 gets merged.
  •  Add more beam problems for documentation.

Week 10: Implementing plotting methods

This week I started working on implementing methods to plot Shear force, bending moment, slope and deflection diagrams. #14967 was created for it.

Mainly four methods were added to the Beam class:

  • plot_shear_force: This method returns a plot for Shear force present in the Beam object.
  • plot_bending_moment: This method returns a plot for Bending moment present in the Beam object.
  • plot_slope: This method returns a plot for slope of the elastic curve of the Beam.
  • plot_delfection: This method returns a plot for the deflection curve of the Beam object.
  • plot_loading_results: This method returns fig object containing subplots of Shear Force, Bending Moment, Slope and Deflection of the Beam object.

Here is a sample notebook demonstrating how to use these plotting methods.

Next Week

  • Make sure #14967 gets merged.
  • Add more beam problems to the documentation.

Week 8 & 9: Beam_3d class

I started implementing Beam_3d class which can be used to find Shear force, Bending moment, Slope, Deflection and other few things for the Beam object.  PR #14883 was created for this.

I implemented Beam_3d class using  this paper as a reference. Actually, like Beam class, it uses a few sets of equations to find certain quantities:

  • To find Shear force and Bending moment
    where [N, Qy, Qz] and [Mx, My, Mz] are the shear force and bending moment along x-y-z-axes respectively (q and m are applied load and moment).
  • To find Slope and Deflection:
    where [wx, wy, wz] and x, θy, θz] are deflection and slope along three axes respectively.

Example for the API:

There is a beam of l meters long. A constant distributed load of magnitude q
is applied along the y-axis from start till the end of the beam. A constant distributed
moment of magnitude m is also applied along the z-axis from start till the end of the beam. Beam is fixed at both of its end. So, deflection of the beam at the both ends
is restricted.

>>> from sympy.physics.continuum_mechanics.beam import Beam_3d
>>> from sympy import symbols
>>> l, E, G, I, A = symbols('l, E, G, I, A')
>>> b = Beam_3d(l, E, G, I, A)
>>> b.apply_support(0, "fixed")
>>> b.apply_support(l, "fixed")
>>> q, m = symbols('q, m')
>>> b.apply_load(q, dir="y")
>>> b.apply_moment_load(m, dir="z")
>>> b.shear_force()
[0, -q*x, 0]
>>> b.bending_moment()
[0, 0, -m*x + q*x**2/2]
>>> b.solve_slope_deflection()
>>> b.slope()
[0, 0, l*x*(-l*q + 3*l*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I)) + 3*m)/(6*E*I)
+ q*x**3/(6*E*I) + x**2*(-l*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I))
- m)/(2*E*I)]
>>> b.deflection()
[0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I))
+ l*m*x**2/(4*E*I) - l*x**3*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I)) - m*x**3/(6*E*I)
+ q*x**4/(24*E*I) + l*x*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I)) - q*x**2/(2*A*G), 0]


As this class is relatively new, it would require a few improvements in the future:

  • As Beam_3d doesn’t use SingularityFunction, I was unable to find a way to represent point load/moments. So for now Beam_3d  only supports continous load (applied over the whole span length of beam).
  • Also, This class assumes that any kind of distributed load/moment is
    applied throughout the span of a beam.

For now, after discussing it with Arihant, we decided to raise NotImplementedError in such cases.

Next Week

  • Make sure PR #14883 gets merge by the end of next week.
  • Start implementing plotting methods for Beam class.

Week 7: Using continuum_mechanics with units

This week I mainly focused on finding and solving a bug due to which continuum_mechanics gave ValueError on using units with the quantities passed. Initially, I created #14856, which included a workaround in the Beam class itself to handle that error. But Arihant suggested opening a separate PR as the error was occurring due to a bug in the separate module.

So, #14865 was created. is_commutative attribute was added in the Prefix class  (setting Prefix.is_commutative to True removed PolynomialError). While doing changes in the PR, another bug appeared:

>>> from sympy.physics.units import meter, newton, kilo
>>> from sympy.physics.units.util import quantity_simplify
>>> quantity_simplify(x*(8*newton + y))
x*(8*newton + y, 1)

This bug was solved with few changes. After #14865 gets merged, continuum_mechanics should work with quantities involving units.

Next Week

  • Make sure #14865 gets merged.
  • Open a Pull Request and start working on 3dbeam class.

Week 6: Max Bending moment and Shear Force

Last week I created #14826 to solve statically indeterminate beam systems using boundary conditions (bc_slope and bc_deflection). This week I mainly focussed on implementing a logic to find the maximum bending moment and maximum shear force in a beam.

Initially, I thought it would be as simple as differentiating the bending_moment and shear force and then solving those using solve. But solve couldn’t represent Interval solutions and hence gave a NonImplemented error, as both of these quantities can occur in Intervals.

So instead of using solve over whole spam length, I found out points of discontinuity in the respective equations using

for term in terms:
    if isinstance(term, Mul):
         term = term.args[-1]

where terms are all the Muls extracted from Add. and term.args[-1] gives us the point of singularity.

Now between two singularity points, our function can be:

  1. A continuous and differentiable function (hence no Interval solution)
  2. or a constant value

for the first scenario, you just use solve over that interval and see values at the endpoint. The higher one of both gives you maxima in that interval. For the second, the constant value is indeed maximum by itself. Then compare maxims of all intervals and return location and its value.

>>> from sympy.physics.continuum_mechanics.beam import Beam
>>> from sympy import symbols
>>> E, I = symbols('E, I')
>>> l, P = symbols('l, P', positive=True)
>>> b = Beam(l, E, I)
>>> R1, R2 = symbols('R1, R2')
>>> b.apply_load(R1, 0, -1)
>>> b.apply_load(R2, l, -1)
>>> b.apply_load(P, 0, 0, end=l)
>>> b.solve_for_reaction_loads(R1, R2)
>>> b.max_bmoment()
(l/2, P*l**2/8)


Next Week

  • Beam class gives ValueError if passed value contains unit. So I would focus on fixing it.
  • Read relevant theory for implementation of 3dBeam class.

Week 5: Solving statically indeterminate beams

This week I worked on solving statically indeterminate beam systems and created #14826 for that. Some work was already done in #14681, which me and Jason reviewed during community bonding period.

Now Beam class uses boundary conditions (bc_deflection and bc_slope) to solve for unknown reactions, hence making statically indeterminate systems solvable.

>>> from sympy.physics.continuum_mechanics.beam import Beam
>>> from sympy import symbols
>>> E, I = symbols('E, I')
>>> F = Symbol('F')
>>> l = Symbol('l', positive=True)
>>> b5 = Beam(l, E, I)
>>> b5.bc_deflection = [(0, 0),(l, 0)]
>>> b5.bc_slope = [(0, 0),(l, 0)]
>>> R1, R2, M1, M2 = symbols('R1, R2, M1, M2')
>>> b5.apply_load(R1, 0, -1)
>>> b5.apply_load(M1, 0, -2)
>>> b5.apply_load(R2, l, -1)
>>> b5.apply_load(M2, l, -2)
>>> b5.apply_load(-F, l/2, -1)
>>> b5.solve_for_reaction_loads(R1, R2, M1, M2)
>>> b5.reaction_loads
{R1: F/2, R2: F/2, M1: -F*l/8, M2: F*l/8}

Next Week

  • Add max_bmoment and mx_shear_force methods to #14826.

Week 4: Solving hinged beams

This week I worked on adding support for beams connected via hinge in #14773. Support for axially fixed beams and its API was implemented last week.

_solve_hinge_beams was added as a helper function to solve such Beams. This method resolves the composite Beam into its sub-beams and then equations of shear force, bending moment, slope and deflection are evaluated for both of them separately. These equations are then solved for unknown reactions and integration constants using the boundary conditions applied on the Beam. Equal deflection of both sub-beams at the hinge joint gives us another equation to solve the system.

So the final API looks like:


>>> from sympy.physics.continuum_mechanics.beam import Beam
>>> from sympy import symbols
>>> l = symbols('l', positive=True)
>>> R1, M1, R2, R3, P = symbols('R1 M1 R2 R3 P')
>>> b1 = Beam(2*l, E, I)
>>> b2 = Beam(2*l, E, I)
>>> b = b1.join(b2,"hinge")
>>> b.apply_load(M1, 0, -2)
>>> b.apply_load(R1, 0, -1)
>>> b.apply_load(R2, l, -1)
>>> b.apply_load(R3, 4*l, -1
>>> b.apply_load(P, 3*l, -1)
>>> b.bc_slope = [(0, 0)]
>>> b.bc_deflection = [(0, 0), (l, 0), (4*l, 0)]
>>> b.solve_for_reaction_loads(M1, R1, R2, R3)
>>> b.reaction_loads
{R3: -P/2, R2: -5*P/4, M1: -P*l/4, R1: 3*P/4}
>>> b.slope().subs(x, 3*l)
>>> b.deflection().subs(x, 2*l)

Next Week

  • See for changes in #14786 to get it merged.
  • Add support for non-horizontal beams.
  • See for any remaining implementation from first two stages of my proposal.