Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added a new method cut_section() which would return a new polygon which lies above the given line that intersects it. #17001

Open
wants to merge 4 commits into
base: master
from

Conversation

Projects
None yet
6 participants
@ishanaj
Copy link
Contributor

commented Jun 9, 2019

A new method cut_section() has been introduced in the Polygon class of the geometry module which returns a segment of the polygon that lies above the given intersecting line.
This would further help in determining the first moment of area of any arbitrary polygon (or a cross-section).

I will be adding tests.

Related PR: #16964

Release Notes

  • geometry
    • Added a method cut_section() for the Polygon class which returns a segment of the polygon that lies above the given intersecting line.

ping @jashan498 @moorepants @oscarbenjamin

Added a new function cut_section() which would return a new polygon (…
…or a segment of the polygon) which lies above the given line that intersects it.
@sympy-bot

This comment has been minimized.

Copy link

commented Jun 9, 2019

Hi, I am the SymPy bot (v147). 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:

  • geometry
    • Added a method cut_section() for the Polygon class which returns a segment of the polygon that lies above the given intersecting line. (#17001 by @ishanaj)

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

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

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

A new method `cut_section()` has been introduced in the `Polygon` class of the `geometry` module which returns a segment of the polygon that lies above the given intersecting line.
This would further help in determining the first moment of area of any arbitrary polygon (or a cross-section).

I will be adding tests.

#### Related PR: #16964 

#### Release Notes

<!-- Write the release notes for this release below. 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 -->
- geometry
    - Added a method `cut_section()` for the Polygon class which returns a segment of the polygon that lies above the given intersecting line.
<!-- END RELEASE NOTES -->

ping @jashan498 @moorepants @oscarbenjamin 

Show resolved Hide resolved sympy/geometry/polygon.py Outdated
improved code to support vertical lines parallel to y axis and also
supporting negetive values of coeff of y
@codecov

This comment has been minimized.

Copy link

commented Jun 9, 2019

Codecov Report

Merging #17001 into master will increase coverage by 11.929%.
The diff coverage is 90%.

@@             Coverage Diff              @@
##           master   #17001        +/-   ##
============================================
+ Coverage   62.22%   74.15%   +11.929%     
============================================
  Files         620      620                
  Lines      160500   160588        +88     
  Branches    37652    37679        +27     
============================================
+ Hits        99864   119076     +19212     
+ Misses      54804    36101     -18703     
+ Partials     5832     5411       -421
@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 9, 2019

The logic used to determine whether a point is above the line is similar to the answer of the math.stackexchange question here

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 9, 2019

I don't understand the purpose of this method.

If you want to find the first moment of area of a polygon is it not better to consider the polygon as being made of triangles?

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 9, 2019

It is kind of related what PR #14434 intends to do.

If you want to find the first moment of area of a polygon is it not better to consider the polygon as being made of triangles?

I haven't considered this method yet. I will have to look into it. Could you please share some source, if possible.
Also, Will the method of triangles will help to calculate the first moment about any axis?

if not intersection_points:
raise ValueError("This line does not intersect the polygon")
points = self.vertices
eq = line.equation()

This comment has been minimized.

Copy link
@smichr

smichr Jun 9, 2019

Member

define x and y as Dummy variables and pass then to equation; there is no point in creating them after the fact below:

x, y = symbols('x y', real=True, cls=Dummy)
eq = line.equation(x, y)
...

This comment has been minimized.

Copy link
@ishanaj

ishanaj Jun 9, 2019

Author Contributor

Just a doubt, what is the significance of _symbols? I used it as it was used in the creation of the line.equation() here

This comment has been minimized.

Copy link
@smichr

smichr Jun 9, 2019

Member

It was a way to handle input from a user for the x and y variables. If they don't provide input then a string is used (e.g. 'x'); if the user provides a symbol then _symbol will not create an error. Since your method doesn't need input from the user, you are free to create your own non-clashing symbol using Dummy.

@@ -783,6 +783,59 @@ def intersection(self, o):
return list(ordered(intersection_result))


def cut_section(self, line):

This comment has been minimized.

Copy link
@smichr

smichr Jun 9, 2019

Member

I think this algorithm will only work for a convex polygon. If there are concave regions you would have to replace all connected points below.

>>> ice = Polygon((0,0),(1,2.5),(2,1),(3,2.5),(4,2),(5,3),(-1,3))
>>> ice.cut_section(Line((0,0),(4.5,3)))
Polygon(Point2D(0, 0), Point2D(9/2, 3), Point2D(1, 5/2), Point2D(3, 5/2), Point2D(-1, 3))

but it should be

>>> ans
Polygon(
 Point2D(0, 0), Point2D(1, 5/2), Point2D(24/13, 16/13), Point2D(12/5, 8/5),
 Point2D(3, 5/2), Point2D(24/7, 16/7), Point2D(9/2, 3), Point2D(-1, 3))

This comment has been minimized.

Copy link
@ishanaj

ishanaj Jun 10, 2019

Author Contributor

@smichr I have made some changes to make it handle concave polygons, please have a look.
The above example gives the correct result now.

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 9, 2019

Any polygon can be represented as a union of triangles. There are algorithms to find the triangles:
https://en.wikipedia.org/wiki/Polygon_triangulation

Simple algorithm: find 3 adjacent vertices A, B, C. These are an "ear" if angle ABC < pi and the line AC does not intersect any other part of the perimeter of the polygon. If so you can separate the triangle ABC from the rest of the polygon by removing vertex B. If the polygon has more than 3 vertices there will always be an ear so you can do this in a loop until the whole polygon is reduced to triangles.

Once you have the polygon represented as triangles all properties such as area, centroid etc can be computed as linear combinations of those of the triangles.

A method to split a polygon into triangles would be of more general use than the cut_section function.

@smichr

This comment has been minimized.

Copy link
Member

commented Jun 9, 2019

consider the polygon as being made of triangles?

And for this there is the code here

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 10, 2019

I will be looking into the Polygon triangulation for determining the first moment of area.
But apart from this, I think it will be a good feature to be added in the geometry module (if it works well).
Would it?


def test_cut_section():
p = Polygon((-1, -1),(1, 2.5),(2, 1),(3, 2.5),(4, 2),(5, 3),(-1, 3))
l = Line((0, 0),(4.5, 3))

This comment has been minimized.

Copy link
@jashan498

jashan498 Jun 11, 2019

Contributor

Using float type for coordinates is causing the tests to fail. You should write them in Rational instead:

p = Polygon((-1, -1),(1, S(5)/2),(2, 1),(3, S(5)/2),(4, 2),(5, 3),(-1, 3))
l = Line((0, 0),(S(9)/2, 3))
assert p.cut_section(l) == Polygon(
Point2D(-9/13, -6/13), Point2D(1, 5/2), Point2D(24/13, 16/13),
Point2D(12/5, 8/5), Point2D(3, 5/2), Point2D(24/7, 16/7),
Point2D(9/2, 3), Point2D(-1, 3), Point2D(-1, -2/3))

This comment has been minimized.

Copy link
@jashan498

jashan498 Jun 11, 2019

Contributor

Similarly here

 Polygon(Point2D(-S(9)/13, -S(6)/13), Point2D(1, S(5)/2),
         Point2D(S(24)/13, S(16)/13), Point2D(S(12)/5, S(8)/5),
         Point2D(3, S(5)/2), Point2D(S(24)/7, S(16)/7),
         Point2D(S(9)/2, 3), Point2D(-1, 3), Point2D(-1, -S(2)/3))
@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 11, 2019

@oscarbenjamin The first moment of area can also be determined by the formula as mentioned in this paper on Pgno. 4 eq-(2.16 and 2.17), which is also used in the current Polygon class to determine the centroid.
Will it be better to use this method or the Polygon triangulation method?

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 11, 2019

There can be lots of ways to calculate the first moment of area. The triangles were just a suggestion because decomposing a polygon into triangles makes it easy to calculate any quantity.

Probably the most interesting question is which methods can be used with symbolic coordinates.

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 11, 2019

I think this method (that I mentioned above) is similar to the current implementations to calculate the centroid and second_moment_of area.
Also, the method of Polygon decomposition seems to be a bit more lengthy, and overwork for calculating just the first moment of area, since polygon decomposition could help us in getting centroid, second moment etc, but we already know those.
So in my humble opinion using this method would be better.
What do you think? @oscarbenjamin

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 11, 2019

That sounds fine.

Actually I didn't realise how many methods were already implemented. Is it not straight-forward to find the first moment of area if you can already compute the centroid and the area?

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 12, 2019

Is it not straight-forward to find the first moment of area if you can already compute the centroid and the area?

Since the formula that I mentioned for the first moment has already been implemented to calculate the centroid, so yes we can directly multiply area and centroid.

But I guess it will give the first moment about the x-axis and y-axis and not about the centroidal (or neutral) axis as we want. Even in current method second_moment_of_area() the formula calculates it about the coordinate axis and then uses parallel axis theorem.

Is there a way similar to parallel axis theorem for first moment?
If not, here's where the cut_section() might help.

@jashan498

This comment has been minimized.

Copy link
Contributor

commented Jun 12, 2019

Shouldn't taking distances y from the centroid of polygon (in the formula ) give first_moment around centroidal axis only?

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 12, 2019

It is better explained in the last page of the notes here, but up to what I had inferred from it and also studied from other sources, y is actually the distance of the centroid of the element dA from the neutral axis. For the sake of simplicity, we can consider dA simply as the part of the polygon above the neutral axis (or any line) and not a small elemental strip.
So we know the location of the neutral axis but we don't know the location of the centroid of the element dA.

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 12, 2019

The parallel axis theorem for first moment of area is straight-forward unless I misunderstand what you mean.

First moment of area around y-axis (line x=0):

Qy = int x dA

First moment around line x=x0:

Qy0 = int (x-x0) dA = Qy - x0*A
@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 12, 2019

Correct me if I am wrong.
I think I just discovered that we don't need the first moment of the whole section, but only the first moment of the part above or below the centroidal axis, since the first moment of the entire area about the centroidal axis will always be zero.
So even if we apply parallel axis theorem, I guess we would get only zero about the centroidal axis

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 12, 2019

I'm not sure what the first moment is used for here apart from finding the centroid. I know that the second moment can come into the calculations:
https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory#Static_beam_equation

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 12, 2019

I think first moment would be used in the calculation of transverse shear stress, if we plan to add stress calculations in beam module.
stress = QV/Ib

@jashan498

This comment has been minimized.

Copy link
Contributor

commented Jun 12, 2019

I think I just discovered that we don't need the first moment of the whole section, but only the first moment of the part above or below the centroidal axis, since the first moment of the entire area about the centroidal axis will always be zero

Actually, this is the way first moment is defined. You take the minimum of int y dA obtained from portions above and below the axis respectively. In case of centroidal axis, both portions will give equal values for int y dA, so you can consider any portion.

You can use this online calculator to verify the results.

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 13, 2019

I guess then we surely need a method to determine the area and the centroid of the section above the neutral axis (or any line) so as to determine the first moment.
IMHO cut_section() is a good way to do so.
I have made a wiki page of how this algorithm works. Please have a look.

Should we be going ahead with the cut_section()?
We can make some changes in it if needed. For eg I think that instead of just returning the polygon segment above, we can also return the polygon segments above as well as below the line. Any suggestions?

@jashan498

This comment has been minimized.

Copy link
Contributor

commented Jun 13, 2019

Should we be going ahead with the cut_section()?

I don’t have any issue withcut_section. It should be helpful in calculating Q at different layers for shear stress.

I think that instead of just returning the polygon segment above, we can also return the polygon segments above as well as below the line.

Yeah, I guess returning a two element list containing both segments will be better than returning just one.

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 13, 2019

Will it also be necessary to implement this method on other shapes e.g. Circle for them to be used in the beam module as well?

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 14, 2019

Will it also be necessary to implement this method on other shapes e.g. Circle for them to be used in the beam module as well?

Yes, I guess a similar functionality would surely be required for ellipses. I will be looking into it

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 14, 2019

For that, we can consider PR #14434
Any suggestions?

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 14, 2019

I think that if we are adding a method like this to Polygon then we should consider whether it can be used by all shapes in the geometry module and how it fits with the general API.

Perhaps another way of thinking about this is as the intersection of the polygon with a half-plane. Using a Line here is strange because the geometry module doesn't distinguish the direction of anti-parallel lines. If we had something like a half-plane represented by a linear inequality like a*x+b*c+d > 0 then we would be able to distinguish which side of the line is the bit we want. That would also generalise better to 3D.

Maybe I'm overcomplicating this. If this particular bit of code just happens to be needed to get the calculations in the beam module working then I think that's fine but it would be a more enticing PR if we could see this method being used for those calculations.

# when coefficient of y is 0, right side of the line is
# considered
compare = (eq.evalf(subs={x: point.x, y: point.y}))/b if b \
else (eq.subs(x, point.x))/a

This comment has been minimized.

Copy link
@oscarbenjamin

oscarbenjamin Jun 14, 2019

Contributor

Can this not just be eq.subs({x: point.x, y:point.y})?

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 14, 2019

Also we should consider how any method added would work with symbolic coordinates. If the concept used is expressed as something well-defined mathematically like an intersection then we have the possibility to use an unevaluated intersection for cases where we don't know which vertices are above or below the line.

Probably in the beam module it would be most useful to get stress/etc as symbolic expressions so you'd want the first moment of area to be a (probably Piecewise) expression in terms of e.g. y.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.