Since you mentioned triangulation, Quarteroni, Sacco, and Saleri's book gives the following formula in the section on multi-dimensional integration: if your domain of integration is a convex polygon $\Omega$ with a triangulation $\mathcal{T}_h$ consisting of $N_T$ triangles, where $h$ is the maximum edge length in $\mathcal{T}_h$, and your quadrature rule has degree of exactness equal to $n$ with all nonnegative weights, then there exists a positive constant $K_n$, independent of $h$, such that the error $E$ is bounded by

$|E| \leq K_n h^{n+1} |\Omega| M_{n+1}$

where $M_{n+1}$ is the maximum value of the modules of the derivatives of order $n+1$ and $|\Omega|$ is the area of $\Omega$. The composite midpoint formula and the composite trapezoidal formula have nonnegative weights and degree of exactness $n=1$ and so we have

$|E| \leq K_1 h^2 |\Omega| M_2$

for some constant $K_1$. It follows that, if your desired error is $\epsilon$, then the maximum length of any edge of your triangulation $h$ must be at most $\epsilon^{1/2}(|\Omega| M_2 K_1)^{-1/2}$, i.e. $O(\epsilon^{1/2}(|\Omega| M_2)^{-1/2})$. If you break $\Omega$ into $N_T$ triangles as I mentioned before, the edge length will be $O(\sqrt{|\Omega|/N_T})$, so you'll need to break $\Omega$ into $O(\frac{|\Omega|^2 M_2}{\epsilon})$ triangles. So, quadratic in the area, linear in the maximum value of the second derivative, and inversely proportional to $\epsilon$. It's not clear to me where convexity enters the picture, so perhaps this will adapt to the general case you've described. The authors state that a proof can be found in Isaacson E. and Keller H. (1966), Analysis of Numerical Methods. Wiley, New York. (This is not my area of expertise so I welcome any corrections)

Since you mentioned triangulation, Quarteroni, Sacco, and Saleri's book gives the following formula in the section on multi-dimensional integration: if your domain of integration is a convex polygon $\Omega$ with a triangulation $\mathcal{T}_h$ consisting of $N_T$ triangles, where $h$ is the maximum edge length in $\mathcal{T}_h$, and your quadrature rule has degree of exactness equal to $n$ with all nonnegative weights, then there exists a positive constant $K_n$, independent of $h$, such that the error $E$ is bounded by

$|E| \leq K_n h^{n+1} |\Omega| M_{n+1}$

where $M_{n+1}$ is the maximum value of the modules of the derivatives of order $n+1$ and $|\Omega|$ is the area of $\Omega$. The composite midpoint formula and the composite trapezoidal formula have nonnegative weights and degree of exactness $n=1$ and so we have

$|E| \leq K_1 h^2 |\Omega| M_2$

for some constant $K_1$. It follows that, if your desired error is $\epsilon$, then the maximum length of any edge of your triangulation $h$ must be at most $\epsilon^{1/2}(|\Omega| M_2 K_1)^{-1/2}$, i.e. $O(\epsilon^{1/2}(|\Omega| M_2)^{-1/2})$. If you break $\Omega$ into $N_T$ triangles as I mentioned before, the edge length will be $O(\sqrt{|\Omega|/N_T})$, so you'll need to break $\Omega$ into $O(\frac{|\Omega|^2 M_2}{\epsilon})$ triangles. So, quadratic in the area, linear in the maximum value of the second derivative, and inversely proportional to $\epsilon$. It's not clear to me where convexity enters the picture, so perhaps this will adapt to the general case you've described. The authors state that a proof can be found in Isaacson E. and Keller H. (1966), Analysis of Numerical Methods. Wiley, New York. (This is not my area of expertise so I welcome any corrections)

UPDATE: I checked the Isaacson book and the relevant section, "7.4 Composite Formulae for Multiple Integrals" doesn't make any mention of convexity, so you may be safe in using the $O(\frac{|\Omega|^2 M_2}{\epsilon})$ that I gave before.