code: def polytope_integrate(
      poly, expr=None, *, clockwise=False, max_degree=None
  ):
      """Integrates polynomials over 2/3-Polytopes.
  
      Explanation
      ===========
  
      This function accepts the polytope in ``poly`` and the function in ``expr``
      (uni/bi/trivariate polynomials are implemented) and returns
      the exact integral of ``expr`` over ``poly``.
  
      Parameters
      ==========
  
      poly : The input Polygon.
  
      expr : The input polynomial.
  
      clockwise : Binary value to sort input points of 2-Polytope clockwise.(Optional)
  
      max_degree : The maximum degree of any monomial of the input polynomial.(Optional)
  
      Examples
      ========
  
      >>> from sympy.abc import x, y
      >>> from sympy import Point, Polygon
      >>> from sympy.integrals.intpoly import polytope_integrate
      >>> polygon = Polygon(Point(0, 0), Point(0, 1), Point(1, 1), Point(1, 0))
      >>> polys = [1, x, y, x*y, x**2*y, x*y**2]
      >>> expr = x*y
      >>> polytope_integrate(polygon, expr)
      1/4
      >>> polytope_integrate(polygon, polys, max_degree=3)
      {1: 1, x: 1/2, y: 1/2, x*y: 1/4, x*y**2: 1/6, x**2*y: 1/6}
      """
      if clockwise:
          if isinstance(poly, Polygon):
              poly = Polygon(*point_sort(poly.vertices), evaluate=False)
          else:
              raise TypeError(
                  "clockwise=True works for only 2-Polytope"
                  "V-representation input"
              )
  
      if isinstance(poly, Polygon):
          # For Vertex Representation(2D case)
          hp_params = hyperplane_parameters(poly)
          facets = poly.sides
      elif len(poly[0]) == 2:
          # For Hyperplane Representation(2D case)
          plen = len(poly)
          if len(poly[0][0]) == 2:
              intersections = [
                  intersection(poly[(i - 1) % plen], poly[i], "plane2D")
                  for i in range(0, plen)
              ]
              hp_params = poly
              lints = len(intersections)
              facets = [
                  Segment2D(
                      intersections[i], intersections[(i + 1) % lints]
                  )
                  for i in range(lints)
              ]
          else:
              raise NotImplementedError(
                  "Integration for H-representation 3D"
                  "case not implemented yet."
              )
      else:
          # For Vertex Representation(3D case)
          vertices = poly[0]
          facets = poly[1:]
          hp_params = hyperplane_parameters(facets, vertices)
  
          if max_degree is None:
              if expr is None:
                  raise TypeError(
                      "Input expression must be a valid SymPy expression"
                  )
              return main_integrate3d(expr, facets, vertices, hp_params)
  
      if max_degree is not None:
          result = {}
          if expr is not None:
              f_expr = []
              for e in expr:
                  _ = decompose(e)
                  if len(_) == 1 and not _.popitem()[0]:
                      f_expr.append(e)
                  elif Poly(e).total_degree() <= max_degree:
                      f_expr.append(e)
              expr = f_expr
  
          if not isinstance(expr, list) and expr is not None:
              raise TypeError(
                  "Input polynomials must be list of expressions"
              )
  
          if len(hp_params[0][0]) == 3:
              result_dict = main_integrate3d(
                  0, facets, vertices, hp_params, max_degree
              )
          else:
              result_dict = main_integrate(
                  0, facets, hp_params, max_degree
              )
  
          if expr is None:
              return result_dict
  
          for poly in expr:
              poly = _sympify(poly)
              if poly not in result:
                  if poly.is_zero:
                      result[S.Zero] = S.Zero
                      continue
                  integral_value = S.Zero
                  monoms = decompose(poly, separate=True)
                  for monom in monoms:
                      monom = nsimplify(monom)
                      coeff, m = strip(monom)
                      integral_value += result_dict[m] * coeff
                  result[poly] = integral_value
          return result
  
      if expr is None:
          raise TypeError(
              "Input expression must be a valid SymPy expression"
          )
  
      return main_integrate(expr, facets, hp_params)
  
imports: from functools import cmp_to_key
  from sympy.abc import x, y, z
  from sympy.core import S, diff, Expr, Symbol
  from sympy.core.sympify import _sympify
  from sympy.geometry import Segment2D, Polygon, Point, Point2D
  from sympy.polys.polytools import LC, gcd_list, degree_list, Poly
  from sympy.simplify.simplify import nsimplify
label_desc: Use image node_id 1 for calling a global function with example usage: polytope_integrate(poly, expr) and returns: main_integrate, result, result_dict, main_integrate3d