Week 7 (30 June – 06 July)

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 @parsoyaarihant 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 (23 June – 29 June)

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 (16 June – 22 June)

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 (9 June – 15 June)

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.

Week 3 (2 June – 8 June)

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
>>> 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.

Week 2 (25 May – 1 June)

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 solvesolveset 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.

Google Summer of Code with SymPy

About a week ago, Google announced the projects selected for this year’s Summer of Code. I am glad to inform that my project, Create a Rich Beam Solving System, for SymPy got selected for GSoC 2018.


SymPy  is a Python library for symbolic mathematics. It aims to become a full-featured Computer Algebra System (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible.


I entered the world of Open Source a few months back and journey since then was full of learning curves, which I enjoyed the most. You learn something new every time you try to contribute. Your code getting reviewed by other contributors, who are usually packed with more experience than you got, is an added advantage. All this inspired me to spend my summer writing code for open source projects and made me apply for GSoC 2018.

Coding period for this year’s GSoC would start on May 14 and should end by August 6. Prior to the coding period is Community Bonding period (April 23 – May 14) which usually is used for learning more about the organization’s community.

I would be updating this blog every week, keeping track of all the progress I make towards the completion of my project.

Community Bonding

The first month of GSoC is the Community Bonding period. In practice, the community bonding period is all about, well, community bonding. Rather than jumping straight into coding, you’ve got some time to learn about your organization’s processes.

During this period I am expected to

  • Discuss the work plan, we would be following this summer, with my mentors and get to know all the tools required for my project.
  • Start a Blog to keep a track of my weekly progress. The Blog is required to have an RSS feed. We are expected to have at least one blog post a week, describing our progress for that week.
  • Sending a pull request with my blog RSS feed to https://github.com/sympy/planet-sympy. For now, I have created this Blog on WordPress but I would be shifting soon to GitHub pages as they look cleaner and show no ads.
  • Review pull requests created by other contributors.

As I am having exams at my university from May 14 until May 26, I would start coding early to cover up for those days. Hopefully, I would be creating my pull request in a day or two.

Till now I have been reviewing #14681, which coincidentally implements one of the points I proposed in my GSoC proposal. Hopefully, I would be creating my first pull request in a day or two.