Using Archimedes-Riemann Theorem to solve higher order integrals 
I am trying to evaluate the integral of $f(x,y) = x + 2y$, where $x,y\in I \equiv [0,2]\times[0,1]$, over all of $I$. The Archimedes-Riemann theorem says that if you have an Archimedean partition of $I$, then you can take the limit of the upper or lower Darboux sums as the number of elements in the partition goes to $\infty$ to evaluate the integral.
Here are some important definition the text gives that I am trying to use:
The Cartesian product of closed intervals is the generalized rectangle, $I$.
If $I = I_1\times I_2 \times \cdots \times I_n$, then for each index $i$ from $1$ to $n$, let $P_i$ be a partition of the $i$th edge $I_i$.
If $J_i$ is an interval in the partition $P_i$, then $J\equiv J_1\times J_2\times \cdots \times J_n$.
$P$ is denoted by $P \equiv (P_1,\dots,P_n)$.
The Lower Darboux Sum is given by $L(f,P_k) = \sum\limits_{J \in P_k} m(f,P_k) vol J$, where $m$ is the minimum of $f(x)$ on $J$.

For my example, $f(x,y) = x+2y, (x,y)\in[0,2]\times[0,1]$.
If I choose $P_k$ to be a partition so that along every edge of $I$, there are $k$ equal segments, then $vol J = 2/k^2$, and the minimum is the value closer to the origin. But I end up getting something like
$L(f,P_k) = \sum\limits_{i=1}^k 6/k^3$ (after some half intelligent hammering away, and that looks like something that will go to $0$. I am just in a Real analysis course and have not taken any topology or Algebra or anything like that.

You seem to be on the right track as you have the right area for each interval in the partition. However, you are not taking into account the value of the function. For each of the intervals in your partition, think about where the minimum value would occur. The functions that you have are increasing. Thus, the minimum value occurs at the bottom left corner. In summation notation
\begin{equation*}
L(f,P_k) = \sum_{i=1}^k \sum_{j=1}^k \frac{2}{k^2} \left[\frac{2}{k}(j-1) + 2\left(\frac{1}{k}(i-1)\right)\right]
\end{equation*}
One way that you can look at the sum is imaging the $i$ running through the $y$ values and the $j$ running through the values of $x$ (I would suggest drawing a picture to see why this is true). Once you distribute the sums and use the well known formula:
$$
1 + 2 + \cdots + n = \frac{n(n+1)}{2}
$$
in the above sums you end up getting 
$$
L(f,P_k) =  \frac{4(k-1)}{k}
$$
Taking a limit (as in the theorem you quote) as $k \to \infty$ we obtain
$$
\int \int_I (x + 2y) dA = \lim_{k \to \infty} L(f,P_k) = 4
$$
If you know anything about calculating double integrals using Fubini's Theorem (http://tutorial.math.lamar.edu/Classes/CalcIII/IteratedIntegrals.aspx) you can see that this is indeed the correct answer.
