Implemented new Method to Calculate Beams with Hinges

[mvg2002]

References to other Issues or PRs

This PR solves the problem described in issue #26570.

Brief description of what is fixed or changed

Implemented a new way to calculate beams with hinges using singularity functions. Using this new method beams with multiple hinges can be calculated.

Other comments

As this is a large pull request where multiple things change at once I will go though it point by point.

1. Added method apply_hinges.
   This method places a hinge anywhere in the beam using a singularity function of order -3.
   It introduces a new unknown: the rotation of the hinge.
   It introduces a new boundry condition: bending moment is 0 at the hinge.

2. Added new property bc_bending_moment which shows the bending moment boundry conditions.
   The setter of this property also adds the boundry conditions to _boundry_conditions.

3. Added the new unknowns and boundry conditions to solve_for_reaction_loads.
   In this method the new unknowns and boundry conditions are added which gives this method the possibility to calculate beams with multiple hinges.

4. Changed the way the join method handles hinges.
   The join method used to handle hinges by setting self._composite_type to hinge. Now that the apply_hinge method is there it can apply a hinge on the new beam at the place where the beams are joined.

5. Changed self._composite_type to self._joined_beam
   As hinges are handled in a different way now it doesn't matter anymore if a beam is joined by a hinge or fixed. For some functions it is still needed to know if a beam is joined or not. That is why self._composite_type is changed to self._joined_beam.

6. Removed _solve_hinge_beams.
   As hinge beams can now be calculated in solve_for_reaction_loads this _solve_hinge_beams method became obsolete.

7. Small change in draw method, so that it can draw multiple hinges if necessary.

8. Added an error message when trying to do I.L.D. calculations on a beam with hinges.
   This was not, and still is not, implemented. This can be done at a later time.

Release Notes

• physics.continuum_mechanics
  • Implemented hinges in a new way which makes it possible to calculate beams with multiple hinges.

[sympy-bot]

✅

Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• physics.continuum_mechanics
  • Implemented hinges in a new way which makes it possible to calculate beams with multiple hinges. (#26594 by @mvg2002)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.14.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
This PR solves the problem described in issue #26570. 


#### Brief description of what is fixed or changed
Implemented a new way to calculate beams with hinges using singularity functions. Using this new method beams with multiple hinges can be calculated. 

#### Other comments
As this is a large pull request where multiple things change at once I will go though it point by point. 

1. Added method `apply_hinges`.
This method places a hinge anywhere in the beam using a singularity function of order -3. 
It introduces a new unknown: the rotation of the hinge. 
It introduces a new boundry condition: bending moment is 0 at the hinge. 

2. Added new property `bc_bending_moment` which shows the bending moment boundry conditions. 
The setter of this property also adds the boundry conditions to `_boundry_conditions`.

3. Added the new unknowns and boundry conditions to `solve_for_reaction_loads`.
In this method the new unknowns and boundry conditions are added which gives this method the possibility to calculate beams with multiple hinges. 

4. Changed the way the `join` method handles hinges. 
The `join` method used to handle hinges by setting `self._composite_type` to hinge. Now that the `apply_hinge` method is there it can apply a hinge on the new beam at the place where the beams are joined. 

5. Changed `self._composite_type` to `self._joined_beam`
As hinges are handled in a different way now it doesn't matter anymore if a beam is joined by a hinge or fixed. For some functions it is still needed to know if a beam is joined or not. That is why `self._composite_type` is changed to `self._joined_beam`. 

6. Removed `_solve_hinge_beams`.
As hinge beams can now be calculated in `solve_for_reaction_loads` this `_solve_hinge_beams` method became obsolete. 

7. Small change in `draw` method, so that it can draw multiple hinges if necessary.

8. Added an error message when trying to do I.L.D. calculations on a beam with hinges. 
This was not, and still is not, implemented. This can be done at a later time. 

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->

* physics.continuum_mechanics
    * Implemented hinges in a new way which makes it possible to calculate beams with multiple hinges. 

<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[mvg2002]

I will also take another look at the doctests that keep failing.

[mvg2002]

I understand that this is a big change. Especially the removal of the _solve_hinge_beams method. I will look out for problems with current functionalities and make sure there are no changes for current users.

For now I will continue with the following 3 points:

1. Change the apply_hinge method in such a way that the rotation jump is also accesible to the user.
2. Add a symbolic test for apply_hinge.
3. Fix the problem in the doctest.

[moorepants]

Especially the removal of the _solve_hinge_beams method.

That's fine. It is a private method so it can be removed without deprecation if user facing functionality doesn't change.

[mvg2002]

After looking at excersice 8 of the beam_problems.rst file, I see that there is a change in output between the new function and the old one. This change happens because the hinge is added to the load equation as a singularity function of order -3. This gives much easier and shorter slope and deflection equations, as I will show below here. But the problem of this change in equations is that the output of for example b.load and b.slope are different than before. That is why the doctest kept failing. I will show the difference below here.

The original output of b.load:

             -2          -1               -1                                -1
      P⋅l⋅<x>     3⋅P⋅<x>     5⋅P⋅<-l + x>                 -1   P⋅<-4⋅l + x>
    - ───────── + ───────── - ────────────── + P⋅<-3⋅l + x>   - ──────────────
          4           4             4                                 2

The new output of b.load:

            2           -3          -2          -1               -1                                -1
      13⋅P⋅l ⋅<-2⋅l + x>     P⋅l⋅<x>     3⋅P⋅<x>     5⋅P⋅<-l + x>                 -1   P⋅<-4⋅l + x>
    - ──────────────────── - ───────── + ───────── - ────────────── + P⋅<-3⋅l + x>   - ──────────────
               48                4           4             4                                 2

It is visible that the load equation became a bit longer because of the addition of the singularity function describing the hinge. The main adventage of this becomes visible when looking at the slope and deflection equation. The slope equation is shown below.

Original output of b.slope:

    ⎛     2               2               2               2⎞               ⎛         1          2               2               2⎞        ⎛         1          2               2               2⎞
    ⎜5⋅P⋅l    P⋅<-2⋅l + x>    P⋅<-3⋅l + x>    P⋅<-4⋅l + x> ⎟           0   ⎜  P⋅l⋅<x>    3⋅P⋅<x>    5⋅P⋅<-l + x>    P⋅<-2⋅l + x> ⎟    0   ⎜  P⋅l⋅<x>    3⋅P⋅<x>    5⋅P⋅<-l + x>    P⋅<-2⋅l + x> ⎟           0
    ⎜────── - ───────────── + ───────────── - ─────────────⎟⋅<-2⋅l + x>    ⎜- ──────── + ──────── - ───────────── + ─────────────⎟⋅<x>    ⎜- ──────── + ──────── - ───────────── + ─────────────⎟⋅<-2⋅l + x>
    ⎝  48           4               2               4      ⎠               ⎝     4          8             8               4      ⎠        ⎝     4          8             8               4      ⎠
    ──────────────────────────────────────────────────────────────────── + ──────────────────────────────────────────────────────────── - ───────────────────────────────────────────────────────────────────
                                    E⋅I                                                                E⋅I                                                                E⋅I

The new output of b.slope:

     ⎛      2           0          1          2               2               2               2⎞
     ⎜13⋅P⋅l ⋅<-2⋅l + x>    P⋅l⋅<x>    3⋅P⋅<x>    5⋅P⋅<-l + x>    P⋅<-3⋅l + x>    P⋅<-4⋅l + x> ⎟    
    -⎜─────────────────── + ───────── - ──────── + ───────────── - ───────────── - ─────────────⎟
     ⎝       48               4          8             8               2               4      ⎠
───────────────────────────────────────────────────────────────────────────────────────────────────
                                                 E⋅I

This is a change that user can experience. If they calculate anything with these functions the result will be the same, but the function as a whole will look different. If a user plots the slope or uses b.slope.subs(x, 3l) the result will be the same, but as shown above the output of b.slope will be different.

My question is if this is something that has to follow depreciation policy. I saw that something has to be depreciated if it's backwards incompatible and user code would stop working after the change. I'm not sure if that is the case here. Knowing all that is stated above here, is it necessary for this code do follow depreciation pollicy?

[moorepants]

My question is if this is something that has to follow depreciation policy.

If the output changes to a mathematically equivalent form, then no deprecation is needed. If the output changes from an incorrect answer to a correct answer, no deprecation is needed. So I think you are fine here.

[mvg2002]

@Tom-van-Woudenberg

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[Tom-van-Woudenberg]

@mvg2002, What happens if you put two hinges at the same place or at the beginning or end of a beam? In the previous implementation that was not possible (or maybe something similar with a beam of length 0).
And what happens if you add too many hinges so that a mechanism emerges? (this was already an option in the old implementation so if this is an issue you might not want to solve it)

[mvg2002]

What happens if you put two hinges at the same place or at the beginning or end of a beam? In the previous implementation that was not possible (or maybe something similar with a beam of length 0).

Putting two hinges at the same place will give a ValueError: duplicate symbols given. This can be handled in a better way. I will add this.

Putting a hinge at the beginning or end of a Beam gives problems, especially when there is a fixed support at the place of the hinge. I will add some error messages for this.

Overall it will probably still be possible to break the code if a user really wants to. It will be quite difficult to catch all possible 'wrong' structures. But you are right that it is better to add the error messages right now if certain problems are known.

And what happens if you add too many hinges so that a mechanism emerges? (this was already an option in the old implementation so if this is an issue you might not want to solve it)

It is still possible to add too many hinges and create a mechanism. If this happens solve_for_reaction_loads will output IndexError: tuple index out of range. This is the output the method gives no matter the problem. This is also described in the following issue comment: #21660 (comment). This issue was quite unrelated to what we are talking about here, so I made a new issue about this #26604.

@Tom-van-Woudenberg

[Tom-van-Woudenberg]

Putting two hinges at the same place will give a ValueError: duplicate symbols given. This can be handled in a better way. I will add this.

Putting a hinge at the beginning or end of a Beam gives problems, especially when there is a fixed support at the place of the hinge. I will add some error messages for this.

Overall it will probably still be possible to break the code if a user really wants to. It will be quite difficult to catch all possible 'wrong' structures. But you are right that it is better to add the error messages right now if certain problems are known.

I think these two possibilities to break the code are relevant to treat with your new implementation as the old method didn't allow you to place multiple hinges at arbitrary points. So I argee with your implementation in 43f3faf

[mvg2002]

@mvg2002 , i think you could think of a better name than new_bending_moment. And could you add an issue for a moment line region of 0?

@Tom-van-Woudenberg, I changed the name to real_bending_moment and made the following issue: #26634

[mvg2002]

@moorepants Regarding the failing checks for symbolic and float inputs. I had never heard about the "garbage in -> garbage out" philosophy, but I understand the idea of it.

Using this idea I think the following checks can be removed:

        if loc in self._applied_rotation_hinges:
            raise ValueError('Cannot place two rotation hinges at the same location.')

        if loc == 0:
            raise ValueError('Cannot place hinge at the beginning of the beam.')

        if loc == self.length:
            raise ValueError('Cannot place hinge at the end of the beam.')

There are also some checks that test less trivial faulty placements:

        if any(loc == arg[1] and arg[2] == -2 for arg in self._applied_loads):
            raise ValueError('Cannot place rotation hinge at the location of a moment load.')

        if any(loc == support[0] and support[1] == 'fixed' for support in self._applied_supports):
            raise ValueError('Cannot place rotation hinge at the location of a fixed support. Change fixed support to pin.')

These lines check for faults a user can make by mistake. This was already possible before the change proposed in PR, but it can be caught now. Do you think these checks can also be removed or is it better to change these checks in a certain way that they can handle float and symbolic inputs?

If you think it is better to keep the tests, do you maybe have some hints on how to do this? I tried a few possible solutions but they did not work.

[Tom-van-Woudenberg]

@Tom-van-Woudenberg, I changed the name to real_bending_moment and made the following issue: #26634

Still a bit vague? What about non_zero_moment?

[mvg2002]

@Tom-van-Woudenberg, I changed the name to real_bending_moment and made the following issue: #26634

Still a bit vague? What about non_zero_moment?

It is not accually the non zero moment. The moment can still be zero at points or even regions in the bending moment plot. The difference is that the bending moment equation no more containes singularity functions of orders lower than 0.

I changed it to non_singular_bending_moment as it is the bending moment equation without the singularity functions that become infinity at one single point. I also added a comment for more clarity.

#Removes the singularity functions of order < 0 from the bending moment equation used in this method
non_singular_bending_moment = sum(arg for arg in self.bending_moment().args if not arg.args[1].args[2] < 0)

[mvg2002]

@Tom-van-Woudenberg and @moorepants, I just pushed a new version without the placement checks when applying a hinge. This is done because of the 'garbage in -> garbage out' phylosophy.

It would be nice if there was a possibilty to place a hinge "just" before or "just" after a support or load. This was not possible in the previous implementation and overshoots the goal of this pull request. I made an issue for this possible improvement, visible here: #26639.

The rest of the given feedback is implemented, are there any more changes necessary?

[moorepants]

I left a spelling suggestion and are you sure it is degrees not radians?

[moorepants]

This is going in! Nice work :)

[mvg2002]

Perfect, thank you!