1) In the syntax (in your file)  both ways of calculating the integral are wrong because you wrote that function  p[t]  depends on  theta  and  r .  In calculating the integral Maple believes that  p[t]  does not depend on  theta  and  r . 

2) If you write explicitly the function  p[t](theta,r)  and  r[5], r[6]  and  theta[tt]  are some constants, it is possible to use both methods. Their efficiency is apparently dependent on the form of the function  p[t](theta,r).

You have an inequality with 3 variables. In general, the solutions will be some spatial regions. Maple  has not any simple description of these regions. But it is easy to find symbolic (or numerical) solution if the two variables are specified.

Example:

Ineq := 3*a[1]*a[2]*a[3]-2*a[1]*a[2]-2*a[1]*a[3]-2*a[2]*a[3]+2*a[1]+2*a[2]+2*a[3]-1 <= 3*a[1]*a[2]*a[3]-((a[1]*a[2]*a[3]+1)/(a[2]*a[3]-a[3]+1)-1)*((a[1]*a[2]*a[3]+1)/(a[1]*a[3]-a[1]+1)-1)-((a[1]*a[2]*a[3]+1)/(a[2]*a[3]-a[3]+1)-1)*((a[1]*a[2]*a[3]+1)/(a[1]*a[2]-a[2]+1)-1)-((a[1]*a[2]*a[3]+1)/(a[1]*a[3]-a[1]+1)-1)*((a[1]*a[2]*a[3]+1)/(a[1]*a[2]-a[2]+1)-1):

solve(eval(Ineq, {a[1] = 3, a[2] = 5}));

simplify(fnormal([evalf(%)]), zero);

eq1 := -3.999168585*10^5*p1*q1^2-2.543619525*10^7*q1*p1^2+7.762255614*10^6*q1^3-3.999168584*10^5*p1^3-0.5000000000e-4*p1+0.1007754033e-3*q1 = 0:

eq2 := -3.999168585*10^5*q1*p1^2-7.762255614*10^6*p1^3+2.543619525*10^7*p1*q1^2-0.5000000000e-4*q1-3.999168584*10^5*q1^3-0.1007754033e-3*p1 = 0:

solve({eq1, eq2});

I did not understand that  pseudo rule means .

bcount:=n->nops(select(t->type(t, odd), [seq(binomial(n,k), k=0..n)])):

Example:

bcount(6);

                             4

It would be interesting to find an explicit formula for this.

Edited.

For any simple polygon (not necessarily convex) area can be found by shoelace formula . The algorithm implemented in Maple in the procedure  Area, which also solves  more general problem - symbolically (or numerically) finds the area of a figure, bounded by a non-selfintersecting piecewise smooth curve. 

convert(sin(x), FormalPowerSeries);

Should be

restart;

Ec := (1/2)*Kc*(2*Pi*h0/Pi/(a0^2+h0^2)-2*Pi*h/Pi/(a^2+h^2))^2;

simplify(Ec);

You have a system of 4 matrix equations in 4 unknowns matrices and matrix parameter y . I think that Maple can not solve the system symbolically, but Maple can solve it numerically for specific values of  B2 , B3  and  y (I took all the square matrix of the third order) . 

I have replaced the names of all the matrices with capital letters and converted the original system in 36 nonlinear scalar equations.

restart;

with(LinearAlgebra):

local D;

B2, B3, Y:=seq(RandomMatrix(3, generator=0..10), i=1..3);

A,B,C,D:=seq(Matrix(3, symbol=i), i=[a,b,c,d]); 

eq2 := A.B+C.D+A:

eq3 := A.C+C.D+C:

eq4 := A.B+C.A+B.C:

eq5 := A.B+A.D+B.C:

fsolve(map(t->Equate(op(t))[], [eq2=B2,eq3=B3,eq4=B2,eq5=Y]));

assign(%):

A,B,C,D;

Expr:=expand((x^2-x-3)^10);

select(t->sign(t)=-1, [op(Expr)]);

nops(%);

restart;

solve({-(1/3)*(eta+5)/(eta-3) = 3*eta/(-2+eta)});

assign(evalf(%[1]));

eta;

Instead of  LinearSolve  should be  LinearAlgebra[LinearSolve] 

If you have  m=1..2  and  x=0..2  then it is a function of 2 variables, so

plot3d(sin(m*x), m = 1 .. 2, x = 0 .. 2, axes=normal);

restart;

taylor(f(x), x = gamma);

subs([x-gamma=e[n], seq((D@@k)(f)(gamma)=k!*c[k]*f(gamma), k=1..5)], %);

 Addition. May be should be  f[k](gamma)  instead of  f(gamma) ?

In addition to the plotting of the ellipse was asked to find a canonical equation. The simplest way to do this - using of the formulas from here

restart;

local D;

f:=(x,y)->3*x^2-3*x*y+6*y^2-6*x+7*y-9;

coeffs(f(x,y));

A,B,C,D,E,F:=%;   # Assigning names to coefficients

theta:=1/2*arctan(B/(A-C));   # Finding the angle of rotation of the canonical system of coordinates relative to the original  coordinate system

solve({-2*A*xc-B*yc=D, -B*xc-2*C*yc=E});   # Finding of the center of the ellipse

assign(%);   # Assigning names to coordinates of the center

# The connection between the canonical  coordinates  and the original  coordinates

x:=xcan*cos(theta)-ycan*sin(theta)+xc;   

y:=xcan*sin(theta)+ycan*cos(theta)+yc;

Eq:=simplify(expand(f(x,y)));   # The equation of the ellipse in canonical coordinates

xcan^2/simplify(sqrt(1/coeff(Eq,xcan^2)*(-tcoeff(Eq))))^`2`+ycan^2/simplify(sqrt(1/coeff(Eq,ycan^2)*(-tcoeff(Eq))))^`2`=1;  # Canonical equation

Visualization:

xcan:=plot(yc+tan(theta)*('x'-xc), 'x'=-2..3.5, color=black):

ycan:=plot(yc-1/tan(theta)*('x'-xc), 'x'=0.1..1.5, color=black):

Ellipse:=plots[implicitplot](f('x','y'),'x'=-2..3.5,'y'=-2..1.5,color=red,thickness=2, gridrefine=5):

labels:=plots[textplot]([[0.4,1.3,"ycan"],[3.2,0.75,"xcan"]], font=[TIMES,ROMAN,14]):

plots[display](xcan,ycan,Ellipse,labels,scaling=constrained);

plot(f(t,a0,b0,c0),  t=0..10);  # or

plot(eval(f(t,a,b,c), {a=a0, b=b0,c =c0), t=0..10);