The double integral of a function f(x,y) over a bounded region
D (of a general shape, not necessarily a rectangle) is defined
simply as the integral over a rectangle R that contains D of
a function F(x,y) that coincides with f(x,y) inside D but
it takes the value 0 outside D.
Since we already know how to integrate over rectangles the above trick is all we need as the definition of the integral over (bounded) regions D of general shape.
Notice, that if f(x,y) is non negative the double integral
over D can still be interpreted as the volume that lies
between D and the surface z=f(x,y). If the boundary of
the region D is not too weird, in particular if
it is of one the forms indicated in the figures below,
the double integrals are given by, 
> R1 := Int(Int(f(x,y),y=g1(x)..g2(x)),x=a..b);
b g2(x) / /   R1 :=   f(x, y) dy dx   / / a g1(x)
this is when the region is like the picture on the left. The other case gives, 
> R2 := Int(Int(f(x,y),x=f1(y)..f2(y)),y=c..d);
d f2(y) / /   R2 :=   f(x, y) dx dy   / / c f1(y)

Compute the double integral of the function, 
> f := (x,y) > x^3*y^2+x*y: 'f(x,y)' = f(x,y);
3 2 f(x, y) = x y + x y
over the region D bounded between the curves, 
> h1 := x^2+x; h2 := x^3x;
2 h1 := x + x 3 h2 := x  x
for non negative x.

We first find the points of intersection of the two curves, h1 and h2. They intersect when x equals: 
> inter := solve(h1h2,x);
inter := 0, 1, 2
The two curves look like this: 
> plot({h1,h2},x=2..2.5,axes=frame);
Since we want only the non negative x's, we need to integrate x from 0 to 2. The integral we need to compute is then, 
> Ans1 := Int(Int(f(x,y),y=h2..h1),x=0..2)=int(int(f(x,y),y=h2..h1),x=0..2);
2 2 x + x / /   3 2 304576 Ans1 :=   x y + x y dy dx =    5005 / / 0 3 x  x
and it evaluates to, 
> ans1 := evalf(304576/5005);
ans1 := 60.85434565

Consider the region D between the circles of radius 1 and 2 in the first quadrant. Look at the picture below: 
> plot({[cos(t),sin(t),t=0..Pi/2],[2*cos(t),2*sin(t),t=0..Pi/2]}, x=3..3,y=3..3);
write the double integral of a general function g(x,y) over the region D
as itereted integrals over x and y.

This region is not of the forms considered above. However, the region D can be decompose into two pieces D1 and D2. The piece D1 is for x between 0 and 1 and D2 when x is between 1 and 2. The integral is then equal to: 
> Ans2 := Int(Int(g(x,y),y=sqrt(1x^2)..sqrt(4x^2)),x=0..1) +
> Int(Int(g(x,y),y=0..sqrt(4x^2)),x=1..2);
2 1/2 2 1/2 1 (4  x ) 2 (4  x ) / / / /     Ans2 :=   g(x, y) dy dx +   g(x, y) dy dx     / / / / 0 2 1/2 1 0 (1  x )