This post is the answer to this question.

The procedure named  IntOverDomain  finds a double integral over an arbitrary domain bounded by a non-selfintersecting piecewise smooth curve. The code of the procedure uses the well-known Green's theorem.

Each section in the border should be specified by a list in the following formats :    
1. If a section is given parametrically, then  [[f(t), g(t)], t=t1..t2]    
2. If several consecutive sections of the border or the entire border is a broken line, then it is sufficient to set vertices of this broken line  [ [x1,y1], [x2,y2], .., [xn,yn] ] (for the entire border should be  [xn,yn]=[x1,y1] ).

Required parameters of the procedure:  f  is an expression in variables  x  and  y , L  is the list of all the sections. The sublists of the list  L  must follow in the positive direction (counterclockwise).

The code of the procedure:

restart;
IntOverDomain := proc(f, L) 
local n, i, j, m, yk, yb, xk, xb, Q, p, P, var;
n:=nops(L);
Q:=int(f,x);  
for i from 1 to n do 
if type(L[i], listlist(algebraic)) then
m:=nops(L[i]);
for j from 1 to m-1 do
yk:=L[i,j+1,2]-L[i,j,2]; yb:=L[i,j,2];
xk:=L[i,j+1,1]-L[i,j,1]; xb:=L[i,j,1];
p[j]:=int(eval(Q*yk,[y=yk*t+yb,x=xk*t+xb]),t=0..1);
od;
P[i]:=add(p[j],j=1..m-1) else
var := lhs(L[i, 2]);
P[i]:=int(eval(Q*diff(L[i,1,2],var),[x=L[i,1,1],y=L[i,1,2]]),L[i,2]) fi;
od; 
add(P[i], i = 1 .. n); 
end proc:

Examples of use.

1. In the first example, we integrate over a quadrilateral:

with(plottools): with(plots):
f:=x^2+y^2:
display(polygon([[0,0],[3,0],[0,3],[1,1]], color="LightBlue"));  # Visualization of the domain of integration
IntOverDomain(x^2+y^2, [[[0,0],[3,0],[0,3],[1,1],[0,0]]]);  # The value of integral

2. In the second example, some sections of the boundary of the domain are curved lines:

display(inequal({{y<=sqrt(x),y>=sin(Pi*x/3)/2,y<=3-x}, {y>=-2*x+3,y>=sqrt(x),y<=3-x}}, x=0..3,y=0..3, color="LightGreen", nolines), plot([[t,sqrt(t),t=0..1],[t,-2*t+3,t=0..1],[t,3-t,t=0..3],[t,sin(Pi*t/3)/2,t=0..3]], color=black, thickness=2));
f:=x^2+y^2: L:=[[[t,sin(Pi*t/3)/2],t=0..3],[[3,0],[0,3],[1,1]], [[t,sqrt(t)],t=1..0]]:
IntOverDomain(f, L);

3. If  f=1  then the procedure returns the area of the domain:

IntOverDomain(1, L);  # The area of the above domain
evalf(%);

IntOverDomain.mw

Edit.