This paper presents numerical solutions to integration problems with bivariate integrands. Using equally spaced nodes in Adaptive Simpson's Rule as a base case, two ways of sampling the domain over which the integration will take place are examined. Drawing from Ouellette and Fiume, Voronoi sampling is used along both axes of integration and the corresponding points are used as nodes in an unequally spaced degree two Newton-Cotes method. Then the domain of integration is triangulated and used in the Triangular Prism Rules discussed by Limaye. Finally, both of these techniques are tested by running simulations over heavily oscillatory and monomial (up to degree five) functions over polygonal regions.