Print[ NIntegrate[ 1/(1 - Cos[x] Cos[y] Cos[z]), {x, 0, Pi}, {y, 0, Pi}, {z, 0, Pi}, Method -> "LocalAdaptive"] / Pi^3 // Timing]
{0.121543, 1.3932012917627028}
Print[ NIntegrate[ Pi Csc[x/2]^2 EllipticK[-Cos[x] Csc[x/2]^4], {x, 0, Pi}] / Pi^3 // Timing]
{0.021917, 1.3932039076998044}
Print[ Gamma[1/4]^4 / (4 Pi^3) // N // Timing]
{0.000029, 1.3932039296856769}
Print[ 4 / Pi^2 EllipticK[1/2]^2 // N // Timing]
{0.000375, 1.3932039296856764}
For `NIntegrate` the default method is `GlobalAdaptive`. You can try different methods and rules to optimize for the function and bounds.
It's a bit of a challenge to get either the first or second forms you provided to converge if requesting 10 or 15 digits of precision. However, your third form involving the Gamma function evaluates to very high numerical precision quickly. I can get 1,000 digits in under 50 milliseconds, which is "good enough"!
Your usage of NIntegrate tuning options shows a few things:
1) That these days I only use Mathematica as a glorified calculator, evoking a mental image of a 500-ton press being used to crack walnuts.
2) Even the best symbolic CAS in the world isn't magic, and requires hand-holding.
3) With the right knowledge even a very hard nut can be easily cracked!
I wonder why Wolfram Research hasn't added some basic multi-threading capabilities to Mathematica to try different "flags" in parallel, racing the various approaches on each CPU core to see which one wins...
It can't even evaluate the answer numerically, it just freezes for minutes!
Edit: about 10 minutes later it spat out 1.393203930, so it can solve the integral, but not easily.
Update: This paper reminds me why I vehemently hate modern mathematical notation, with subscripts and superscripts can mean anything at all and everyone's idea of "convention" is different. On the second page there's an algebraic expression for the solution, which at first glance seems to be:
4/Pi^2 * EllipticK[Sqrt[1/2 - Sqrt[1 - w^(-2)]/2]]^2
But this evaluates to 1.76351 not 1.39320 as given in the paper. The "K^2" in the paper is not the "complete elliptic integral of the first kind" squared, it's some other function that squares its input somehow. Or something. But that begs the question of why formulate this in terms of a function of a square root squared!?
Helpfully, they provide the definition in terms of the hypergeometric function a bit lower down, but that appears to be wrong as well, providing yet another result if evaluated numerically.