Computing an integral in $\mathbb{R}^2$ by definition Given the function $f(x,y) = xy$, how do I calculate $$\int_{[0,1]\times [0,1]}f(x,y)$$ by the definition of the Riemann integral, without showing it is Riemann integrable, since this comes from continuity on f?
Here are the definitions to be used:
Jordan content: If $ I = [a_1,b_1]\times[a_2,b_2]\times\ldots\times[a_N,b_N]\subset\mathbb{R}^N $ is a compact interval, the its Jordan content, $\mu(I)$, is $$ \mu(I) = \prod_{j=1}^N(b_j-a_j)  $$
Riemann Sum: Let $I\subset\mathbb{R}^N$ be a compact interval. Let $f:I\rightarrow\mathbb{R}^K$ be a function. Let $P$ be a partition of $I$ into $I_\nu$. For each $\nu$, choose some $x_\nu\in I_\nu$. Then the Riemann sum of $f$ corresponding to $P$ is $$  S(f,P) = \sum_{\nu}f(x_\nu)\mu(I_\nu) $$
Riemann Integral: Let $I\subset\mathbb{R}^N$ be a compact interval. Let $f:I\rightarrow\mathbb{R}^K$ be a function. Suppose that there is $y\in \mathbb{R}^N$ such that for each $\epsilon>0$, there is a partition $P_\epsilon$ of $I$ such that for each refinement $P$ of $P_\epsilon$ and for any Riemann sum $S(f,P)$ corresponding to $P$, we have $\|S(f,P)-y\|<\epsilon$. Then $f$ is said to be Riemann integrable on $I$, and $y$ is called the Riemann integral of $f$ over $I$, so $$y = \int_{I}f$$
To calculate the integral, do I need to find a specific partition of $[0,1]\times[0,1]$ into $I_\nu$, choose a specific $x_\nu$ for each $I_\nu$, then compute the Riemann sum? If so, how do I know I'm choosing a good partition? 
 A: I will use the Darboux integral, which is equivalent to the Riemann integral:
Let $\Omega=[0,1]\times [0,1].$ First, we'll construct a partition of $\Omega$:
Let $n,m\in\mathbb{N}$. Then, we define a partition of $\Omega$ by $$P_{n,m}=\{x_0,x_1,\cdots,x_n\}\times\{y_0,y_1,\cdots,y_m\}=\left\lbrace0,\frac{1}{n},\frac{2}{n},\cdots,\frac{n-1}{n},1\right\rbrace\times\left\lbrace0,\frac{1}{m},\frac{2}{m},\cdots,\frac{m-1}{m},1\right\rbrace.$$ Let $i=1,\cdots,n$ and $j=1,\cdots,m$. Define $R_{i,j}=[x_{i-1},x_i]\times[y_{j-1},y_j]$, with $\Omega=\bigcup\limits_{i,j} R_{i,j}.$ Note that $$M_{i,j}(f)=\sup\limits_{(x,y)\in R_{i,j}}|f(x,y)|=x_iy_j,$$ and $$m_{i,j}(f)=\inf\limits_{(x,y)\in R_{i,j}}|f(x,y)|=x_{i-1}y_{j-1}.$$ Note further that content of each $R_{i,j}$ is $\Delta_{i,j}=\frac{1}{nm}.$ Now, we can calculate the upper and lower sums as $$U(f,P_{n,m})=\sum\limits_{i=1}^n\sum\limits_{j=1}^m M_{i,j}(f)\Delta_{i,j}=\sum\limits_{i=1}^n\sum\limits_{j=1}^m \frac{i}{n}\frac{j}{m}\frac{1}{nm}=\frac{1}{4}+\frac{1}{4m}+\frac{1}{4n}+\frac{1}{4nm}$$ and
 $$L(f,P_{n,m})=\sum\limits_{i=1}^n\sum\limits_{j=1}^m m_{i,j}(f)\Delta_{i,j}=\sum\limits_{i=1}^n\sum\limits_{j=1}^m \frac{(i-1)}{n}\frac{(j-1)}{m}\frac{1}{nm}=\frac{1}{4}-\frac{1}{4m}-\frac{1}{4n}+\frac{1}{4nm}.$$ Now, let $\epsilon>0$ and  choose $N,M\in\mathbb{N}$ so that $N,M>\frac{1}{\epsilon}$, and define the partition $P$ as $P_{N,M}$, which is given in the manner constructed earlier. Then,
 $$U(f,P)-L(f,P)=\frac{1}{2N}+\frac{1}{2M}<\frac{\epsilon}{2}+\frac{\epsilon}{2}=\epsilon. $$ Thus, $f$ is Darboux integrable on $\Omega$.
 Further, we have that the infimum of the upper sum and the supremum of the lower sum both occur as $n,m\rightarrow\infty$; hence, they converge to $1/4$, which means that this is the value of the integral. If we want to be explicit, $$L(f,P)\leq\int\limits_\Omega f(x,y)\ d\Omega\leq U(f,P),$$ which implies that $$\frac{1}{4}-\frac{\epsilon}{4}-\frac{\epsilon^2}{4}\leq\int\limits_\Omega xy\ d\Omega\leq \frac{1}{4}-\frac{\epsilon}{4}+\frac{\epsilon^2}{4}.$$ Since $\epsilon>0$ was arbitrary, we conclude that $$\frac{1}{4}\leq\int\limits_\Omega xy\ d\Omega\leq \frac{1}{4},$$ from which it follows that $$\int\limits_\Omega xy\ d\Omega=\frac{1}{4}.$$
