How can I reduce the error estimate of numerical integration in Mathematica?

Question

Before introducing the integral I want to go through some definitions.

Define triangles $T_1$ and $T_2$ by the points $\{a,b,c\}$ and $\{d,e,f\}$ respectively.

Define $$D(u,v,x,y) = sign(\det (y-u,y-v,y-x)),$$

and define $$n_\cap (T_1, x, y) = \frac{1}{8} (D(a, b, x, y) + D(b, c, x, y) + D(c, a, x, y) + \\D(a, b, x, y)D(b, c, x, y)D(c, a, x, y)) (1 - D(a, b, c, d)D(a, b, c, e)).$$

Say $lk(T_1,T_2)^2$ is defined as

$$lk(T_1, T_2)^2 = [n_\cap (T_1, d, e) + n_\cap(T_1, e, f) + n_\cap (T_1, f, d)]^2.$$

How can I reduce the error of numerical integration of the following integral? $$\int_{\Omega^3} \int_{\Omega^3} lk(T_1, T_2)^2,$$

where $\Omega = [0,1]^3$ is the space in which we take the points $a,b,c,d,e,f$ to generate the triangles $T_1$ and $T_2$. We integrate over the points $\{a, b, c\}$ and $\{d, e, f\}$ (hence the $\Omega^3$). Note that the points $a,b,c,d,e,f$ are points in the unit box.

I have tried to calculate the value of this integral numerically using Mathematica and the result is approximately 0.15... the error estimate is around 0.0016042 and I haven't been able to reduce this error.

I am not used to work with Mathematica but I have tried to use Global and Local adaptive strategies and haven't been succesful. I have also tried setting an Acuracy Goal, Max and Min recursions but this has not worked at all. I even tried to change the integration method but this has not worked either.

Any advice of how to reduce the error of the numerical integration or how to calculate the integral symbolically would be appreciated.

The code for the integral in Mathematica is

a = {a1, a2, a3};
b = {b1, b2, b3};
c = {c1, c2, c3};
d = {d1, d2, d3};
e = {e1, e2, e3};
f = {f1, f2, f3};
x = {x1, x2, x3};
y = {y1, y2, y3};

lk2 := ((1/
       8 (Sign[Det[{e - a, e - b, e - d}]] + 
        Sign[Det[{e - b, e - c, e - d}]] + 
        Sign[Det[{e - c, e - a, e - d}]] + (Sign[
           Det[{e - a, e - b, e - d}]]*
          Sign[Det[{e - b, e - c, e - d}]]*
          Sign[Det[{e - c, e - a, e - d}]])) (1 - (Sign[
           Det[{d - a, d - b, d - c}]]*
          Sign[Det[{e - a, e - b, e - c}]]))) + (1/
       8 (Sign[Det[{f - a, f - b, f - e}]] + 
        Sign[Det[{f - b, f - c, f - e}]] + 
        Sign[Det[{f - c, f - a, f - e}]] + (Sign[
           Det[{f - a, f - b, f - e}]]*
          Sign[Det[{f - b, f - c, f - e}]]*
          Sign[Det[{f - c, f - a, f - e}]])) (1 - (Sign[
           Det[{e - a, e - b, e - c}]]*
          Sign[Det[{f - a, f - b, f - c}]]))) + (1/
       8 (Sign[Det[{d - a, d - b, d - f}]] + 
        Sign[Det[{d - b, d - c, d - f}]] + 
        Sign[Det[{d - c, d - a, d - f}]] + (Sign[
           Det[{d - a, d - b, d - f}]]*
          Sign[Det[{d - b, d - c, d - f}]]*
          Sign[Det[{d - c, d - a, d - f}]])) (1 - (Sign[
           Det[{f - a, f - b, f - c}]]*
          Sign[Det[{d - a, d - b, d - c}]]))))^2

NIntegrate[lk2, {a1, 0, 1}, {a2, 0, 1}, {a3, 0, 1}, {b1, 0, 1}, {b2, 
  0, 1}, {b3, 0, 1}, {c1, 0, 1}, {c2, 0, 1}, {c3, 0, 1}, {d1, 0, 
  1}, {d2, 0, 1}, {d3, 0, 1}, {e1, 0, 1}, {e2, 0, 1}, {e3, 0, 1}, {f1,
   0, 1}, {f2, 0, 1}, {f3, 0, 1}]
```

Answer 1

On such a high dimensional integral, the default rule is the Monte Carlo rule. You can increase the number of points. I also increased the PrecisionGoal, so that the error estimate will be reported.

NIntegrate[..., 
 Method -> {"MonteCarloRule", "Points" -> 10^6}, PrecisionGoal -> 6]

NIntegrate::maxp: The integral failed to converge after 3000000 integrand evaluations. NIntegrate obtained 0.15226900000000002 and 0.0002590579838772007 for the integral and error estimates.

(*  0.152269  *)

The Monte Carlo Rule error is proportional to $1/\sqrt{N}$ where $N$ is the number of "Points" (under certain assumptions). It converges slowly. To get another digit of precision, use about 100 times as many points and wait about 100 times as along.