I need to compute a certain integral $$\int_{-h}^0\int_{-h}^0 f(x,y).$$ I've split the square $[-h,0] \times [-h,0]$ into $N \times N$ equally spaced subsquares. Each subsquare has been split into two triangles along the main diagonal and for each triangle passing through points $$(\theta_i, \theta_j),\quad (\theta_i, \theta_j + h),\quad (\theta_i + h, \theta_j)$$ function $f(x,y)$ has been replaced with the plane passing through 3 points $$\big((\theta_i, \theta_j), f(\theta_i,\theta_j)\big), \big((\theta_i + h, \theta_j),\quad f(\theta_i + h,\theta_j)\big), \big((\theta_i, \theta_j + h),\quad f(\theta_i,\theta_j + h)\big).$$
Now I'm trying to find a good estimate in error after replacing original integral with the so called "piecewise-linear-on-each-triangle". (I don't know how my procedure is called). I thought that it reminds me of the integration with the trapezoidal rule because they coincide for 1-dimensional integrals.
I'm trying to find any information about how a procedure for approximating double integrals would look like if it has been generalized from the trapezoidal rule integration for one-dimensional intergrals (https://en.wikipedia.org/wiki/Trapezoidal_rule#CITEREFAtkinson1989) and I've tried to generalize the proof for error estimate given on wiki but it did not work out for double integrals.
Especially what I'm interested about is the error estimation and a proof behind it. Bonus points if one can say something about my procedure specifically.
Do you specifically want to use the triangles and the piecewise linear interpolator? You can easily generalize the trapezoidal rule for double integrals just ny applying it to the iterated integrals: \begin{align*} \int_{-h}^0 \int_{-h}^0 f(x,y) dx dy \approx & \int_{-h}^0 \sum_{j=0}^N w_j f(\theta_j ,y) dy=\sum_{j=0}^N w_j \int_{-h}^0 f(\theta_j,y)dy\\ \approx & \sum_{j=0}^Nw_j\sum_{i=0}^N w_i f(\theta_j,\theta_i) = \sum_{i,j=0}^N W_{ij}f(\theta_j, \theta_i) \end{align*}
where $w_0= w_N = \frac{h}{2}$, $w_i = h, i = 1, \cdots,N-1$ and $W_{ij} = w_i w_j$. In this situation you'll find that the error formula comes directly from the 1d case.
As for you specific situation, you may want to search for the "finite element method", that heavily relies on numerical integration on partitions of the domain (namely triangles) to discretize differential equations.