How to calculate the total area of $frac{cos (2 y) cos (2 x-y)+cos (3 y)}{cos (2 x-y)+cos (y)}$ enclosed in the range $[0,1]$?

You calculate the length of a curve $f(x)$ between $a$ and $b$ by $int_a^b sqrt{1+f'(x)^2}dx$, see https://en.wikipedia.org/wiki/Arc_length.

The 2d generalisation would logically be $int_A sqrt{1+(partial_x f(x,y))^2+(partial_y f(x,y))^2}dA$. Implementing in mathematica would go like

f[x_,y_] := (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + Cos[y])    
R = ImplicitRegion[{0 < x < 2 Pi, 0 < y < 2 Pi, 0 < f[x, y] < 1}, {x, y}]
NIntegrate[Sqrt[1 + D[f[x, y], x]^2 + D[f[x, y], y]^2], {x, y} [Element] R]

This gives an answer of about 9062.5, which seems to me like like too high of a number, so I might have made a mistake.

Edit: Second try. This answer seems more reasonable and agrees with the other solutions. Some error messages are still given about convergence about which I know very little.

integrand = Sqrt[1 + D[f[x, y], x]^2 + D[f[x, y], y]^2]
NIntegrate[If[0 < f[x, y] < 1, integrand, 0], {x, 0, 2 Pi}, {y, 0, 2 Pi}]

Edit: I thought you were interested in the area of the surface of your function. If you were interested in the are of the domain where $0<f(x,y)<1$, I misunderstood, and you should follow Daniel's answer.

Let us use a straightforward approach.

f = (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + Cos[y])

Its plot is built by

Plot3D[f,{x,0,2*Pi},{y,0,2*Pi}]

The result of

r = ImplicitRegion[0 < (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + 
Cos[y]) && (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + 
Cos[y]) < 1 && x >= 0 && x <= 2*Pi && y >= 0 && y <= 2*Pi, {x, y}]
Region[r]

shows the region of the integration. Next,

FunctionPeriod[f, x]

Pi

and

FunctionPeriod[f,y]

2*Pi

Hence, the required area approximately equlals

2*NIntegrate[Sqrt[1 + D[f, x]^2 + D[f, y]^2]*
Boole[0 < (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + 
Cos[y]) && (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] +
Cos[y]) < 1], {x, 0, Pi}, {y, 0, 2*Pi}, PrecisionGoal -> 2, 
AccuracyGoal -> 2, WorkingPrecision -> 15, Method -> "GlobalAdaptive", MaxRecursion -> 30]

39.6109795138184

and a warning "NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small."

Following Lukas formula for the area, you can use NIntegrate with HeavisideTheta functions to define the domain. You can add accuracy and/or precision options if required.

  NIntegrate[
     Sqrt[1 + D[f[x, y], x]^2 + D[f[x, y], y]^2] HeavisideTheta[
       f[x, y]] HeavisideTheta[1 - f[x, y]], {x, 0, 2 Pi}, {y, 0, 2 Pi}]
    (*39.33*)

Too long for a comment:

Notice: Hacky hack. I like to see what the function looks like, so...

Plot3D[ (Cos[2 y] Cos[2 x - y] + Cos[3 y])/( Cos[2 x - y] + Cos[y]), {x, 0, 2 π}, {y, 0, 2 π}, PlotRange -> {0, 1}, ClippingStyle -> None, PlotPoints -> {100, 100}, ColorFunction -> Function[{x, y, z}, Hue[z]]]

So, it's pretty clear the surface area is more than 6 (2 pi). Either use Lucas's formula, or define the region and use the surface area command built into Mathematica. For the hackiest hack, just double the surface area of the second and third objects from the top