Repository: sympy/sympy

Issue:
`polytope_integrate`: Polygon vertex order assumption
The `Polygon` class has a positive area if the vertices are in ccw order. In contrast the default behaviour of `polytope_integrate` is given the same polygon that the integral of, for example, the one function is negative.

This inconsistent behaviour is especially noticeable because the `polytope_integrate` docs refer to the `Polygon` class docs which state "ccw => positive area". Only the API section gives an implicit statement that cw order is assumed.

In addition the reordering algorithm is apparently broken on debian 9 + python3.5.3 + sympy1.3
```
Traceback (most recent call last):
  File "/home/xxx/yyy/jump2d.py", line 41, in <module>
    B = simplify(polytope_integrate(_D, p0, clockwise=True))
  File "/home/xxx/.local/lib/python3.5/site-packages/sympy/integrals/intpoly.py", line 67, in polytope_integrate
    poly = point_sort(poly)
  File "/home/xxx/.local/lib/python3.5/site-packages/sympy/integrals/intpoly.py", line 1012, in point_sort
    n = len(poly)
TypeError: object of type 'Polygon' has no len()
```



Discussion:

Patch:
diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py
index 7aa41207ac..776d2d87f8 100644
--- a/sympy/integrals/intpoly.py
+++ b/sympy/integrals/intpoly.py
@@ -62,9 +62,9 @@ def polytope_integrate(poly, expr=None, **kwargs):
     clockwise = kwargs.get('clockwise', False)
     max_degree = kwargs.get('max_degree', None)
 
-    if clockwise is True:
+    if clockwise:
         if isinstance(poly, Polygon):
-            poly = point_sort(poly)
+            poly = Polygon(*point_sort(poly.vertices), evaluate=False)
         else:
             raise TypeError("clockwise=True works for only 2-Polytope"
                             "V-representation input")
@@ -1009,19 +1009,20 @@ def point_sort(poly, normal=None, clockwise=True):
     >>> point_sort([Point(0, 0), Point(1, 0), Point(1, 1)])
     [Point2D(1, 1), Point2D(1, 0), Point2D(0, 0)]
     """
-    n = len(poly)
+    pts = poly.vertices if isinstance(poly, Polygon) else poly
+    n = len(pts)
     if n < 2:
-        return poly
+        return list(pts)
 
     order = S(1) if clockwise else S(-1)
-    dim = len(poly[0])
+    dim = len(pts[0])
     if dim == 2:
-        center = Point(sum(map(lambda vertex: vertex.x, poly)) / n,
-                       sum(map(lambda vertex: vertex.y, poly)) / n)
+        center = Point(sum(map(lambda vertex: vertex.x, pts)) / n,
+                        sum(map(lambda vertex: vertex.y, pts)) / n)
     else:
-        center = Point(sum(map(lambda vertex: vertex.x, poly)) / n,
-                       sum(map(lambda vertex: vertex.y, poly)) / n,
-                       sum(map(lambda vertex: vertex.z, poly)) / n)
+        center = Point(sum(map(lambda vertex: vertex.x, pts)) / n,
+                        sum(map(lambda vertex: vertex.y, pts)) / n,
+                        sum(map(lambda vertex: vertex.z, pts)) / n)
 
     def compare(a, b):
         if a.x - center.x >= S.Zero and b.x - center.x < S.Zero:
@@ -1054,9 +1055,7 @@ def compare3d(a, b):
         elif dot_product > S.Zero:
             return order
 
-    if dim == 2:
-        return sorted(poly, key=cmp_to_key(compare))
-    return sorted(poly, key=cmp_to_key(compare3d))
+    return sorted(pts, key=cmp_to_key(compare if dim==2 else compare3d))
 
 
 def norm(point):


Test patch:
diff --git a/sympy/integrals/tests/test_intpoly.py b/sympy/integrals/tests/test_intpoly.py
index 01b00f171c..02c0611313 100644
--- a/sympy/integrals/tests/test_intpoly.py
+++ b/sympy/integrals/tests/test_intpoly.py
@@ -3,11 +3,11 @@
 from sympy.core import S
 
 from sympy.integrals.intpoly import (decompose, best_origin,
-                                     polytope_integrate)
+                                     polytope_integrate, point_sort)
 
 from sympy.geometry.line import Segment2D
 from sympy.geometry.polygon import Polygon
-from sympy.geometry.point import Point
+from sympy.geometry.point import Point, Point2D
 from sympy.abc import x, y, z
 
 
@@ -490,6 +490,14 @@ def test_polytope_integrate():
          y * z: 3125 / S(4), z ** 2: 3125 / S(3), y ** 2: 3125 / S(3),
          z: 625 / S(2), x * y: 3125 / S(4), x ** 2: 3125 / S(3)}
 
+def test_point_sort():
+    assert point_sort([Point(0, 0), Point(1, 0), Point(1, 1)]) == \
+        [Point2D(1, 1), Point2D(1, 0), Point2D(0, 0)]
+
+    fig6 = Polygon((0, 0), (1, 0), (1, 1))
+    assert polytope_integrate(fig6, x*y) == S(-1)/8
+    assert polytope_integrate(fig6, x*y, clockwise = True) == S(1)/8
+
 
 def test_polytopes_intersecting_sides():
     fig5 = Polygon(Point(-4.165, -0.832), Point(-3.668, 1.568),