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.

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.

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, Q_y, Q_z] and [M_x, M_y, M_z] 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 [w_x, w_y, w_z] 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.

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.

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.

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.

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)
-7*P*l**2/(48*E*I)
>>> b.deflection().subs(x, 2*l)
7*P*l**3/(24*E*I)

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.

This week me, Arihant and Jason discussed the API for creating composite Beam objects. Initially, we used Piecewise and SingularityFunction to represent our changing second_moment but then we agreed upon .join method to represent such beams. So the final API was like:

>>> b1 = Beam(2, E, 1.5*I)
>>> b2 = Beam(2, E, I)
>>> b = b1.join(b2, "fixed")
>>> b.length
4
>>> b.second_moment
Piecewise((1.5*I, x <= 2), (I, x <= 4))

Here b1.join(b2, "fixed") joins b2 at the right end of b1 via a fixed connection.All this was implemented in #14773 and hopefully it would be merged in coming few days.

I also created #14786 at the end of this week implementing apply_support and max_deflection methods.

apply_support is an easier way to apply support structures on our Beam object rather than adding all the reaction loads and moments and constraints on it by yourself. Its API is not finalised yet but for now it is something like:

>>>b.apply_support(position, type='hinge')

where position represents the position at which support was applied

and type is type of support structure. It can be either hinge, roller or cantilever.

Next Week

• Add support for composite beams connected via hinge.
• Add support for non-horizontal beams.
• See for any remaining implementation from first two stages of my proposal.

Exams kept me occupied for previous 11 days so I wasn’t able to contribute much. Coming to this week’s work, I created two pull request:

• #14751 :This PR implemented remove_load method to remove previously applied loads on the beam object. This method is  little different from adding a negative load to make net equal to zero and would work only if that particular load exists on beam.
•  #14753  : This PR implemented point_cflexure method to find point of contraflexure.

#14751 also included applied_loads method which keeps a track of all load applied on beam.

It is different from b.load as it treat each load as a separate entity. load property would sum up all the loads at a particular point but applied_loads will still show them as separate loads. For example

>>> b.apply_load(4, 2, -1)
>>> b.apply_load(2, 2, -1)
>>> b.load
6*SingularityFunction(x, 2, -1)
>>> b.applied_loads
[(4, 2, -1, None), (2, 2, -1, None)]

The only difficulty with #14753  occured in finding solution of moment_curve .  Actually moment is zero outside the spam length too, which made solve to return a solution in form of Interval which is not a compatible return type for solve. solveset was of no help too as it can’t be used with multivariate expressions. The problem, however, was solved by wrapping bending_moment  with a Piecewise function with its value equal to float("nan") outside the spam length.

Next Week

• Make beam.py compatible to solve non-prismatic beams.
• Try to find a way around for the issue occurring in PR #14681.
• Support for non-horizontal beams.