It's a very difficult task to compute u3 even at a single point.
For example, u3(0) can be simplified to the triple integral of
f := -(2.533317577*10^7*I)*z*exp((-.64-12.56637062*I)*x^2-(7.853982892*10^7*I)*z^2-(25.13274124*I)*y^2)*y*x*BesselJ(0, 25.13274123*y*x)*BesselJ(0, 62831.85308*z*y)

But f is very oscillatory and to compute its triple integral over [0,1]^3 is almost impossible even using MonteCarlo method.
Probably some special methods will be needed.
Just look at this plot:
plot( eval(Re(f),[x=1/2,y=1/2]),z=0..1 );

P.S. Was this integral invented, or it appears in a concrete problem?

eq:=proc(L::listlist,x::name,y::name,z::name)
  primpart(LinearAlgebra:-Determinant(<Matrix([[x,y,z],L[]])|<1,1,1,1>>))=0
end:

map(eq,L,x,y,z);

TriOnUnitSphere:=proc(A::Vector[column](3),B::Vector[column](3),C::Vector[column](3)) #A,B,C = directions
uses LinearAlgebra,plots,plottools;
local a:=Normalize(A,2),b:=Normalize(B,2),c:=Normalize(C,2),u,v;
display(sphere(),  plot3d(1.01*Normalize(c+v*(a+u*(b-a)-c),2), u=0..1,v=0..1,color=red,style=surface,_rest),
pointplot3d([1.01*a,1.01*b,1.01*c],symbolsize=20,color=blue,symbol=solidsphere)  )
end:

TriOnUnitSphere(<1,2,3>, <-1,1,-2>, <-1,-1,3>);

TriOnUnitSphere(<1,0,0>, <-10,1,0>, <0,1,5>, grid=[200,200]);

with(LinearAlgebra):with(plots):with(plottools):
A1:=<0,0,1>: A2:=<0,3/5,4/5>: A3:=<4/5,3/5,0>: A4:=<5/13,0,12/13>:
seg:= (a,b)->spacecurve(Normalize(a+t*(b-a),2),t=0..1, thickness=3, color=red):
display(sphere([0,0,0],1),seg(A1,A2),seg(A2,A3),seg(A3,A4),seg(A1,A4));

You can use e.g.

binaryvariables = indets(L, specindex(integer,x));

where L is a list containing the constraints (and the function).
Or, you may generate their sequence using seq.
 

f:=(x,y) ->exp(-x^2-y^2):
L:=4: n:=64:
r1:=rand(0..1.0): r:=rand(-1..1.0):
a:=Vector(n,L*r): b:=Vector(n,L*r):
c:=Vector(n,0.3*r):
da:=Vector(n, i->0.5*''r1()''):
P:=proc(aa)
  global n,a,b,c,da,h,L;
  a:=a+eval~(da);
  map[inplace](frem,a,L);
  h:=add(c[i]*f(x+a[i], y+b[i]),i=1..n):
  plot3d(h, x=-L..L,y=-L..L, color="Aquamarine", scaling=constrained, view=-2..2, axes=none, style=surface, lightmodel=light1);
end:
plots:-display(seq(P(aa),aa=1..50),insequence=true);
## or ...
Explore(P(aa), aa=0..1, animate,loop,size=[800,800]);

Download req-int.mw

restart;
f:=r^2 *cos(theta)+r*sin(theta):
r1:= 0.3+0.1*cos(theta):
r2:= 0.5+0.1*cos(theta):
plots:-densityplot(f, theta=0..2*Pi, r=r1 .. r2,  
coords=polar, colorstyle = HUE, style = patchnogrid);

If you need labels:

plots:-display(%, labels=["x","y"]);

c:=Matrix([    # result = 22.613, [4, 2, 5, 3, 6, 1, 4]
[0., 3.60555127546399, 2.00000000000000, 8.06225774829855, 5.09901951359278, 1.41421356237310],
[3.60555127546399, 0., 2.23606797749979, 6.32455532033676, 2.23606797749979, 3.60555127546399],
[2.00000000000000, 2.23606797749979, 0., 8.06225774829855, 3.16227766016838, 1.41421356237310],
[8.06225774829855, 6.32455532033676, 8.06225774829855, 0., 8.06225774829855, 9.00000000000000],
[5.09901951359278, 2.23606797749979, 3.16227766016838, 8.06225774829855, 0., 4.47213595499958],
[1.41421356237310, 3.60555127546399, 1.41421356237310, 9.00000000000000, 4.47213595499958, 0.]]):

n := sqrt(numelems(c)):  N:={seq(1..n)}:
z:=add(add( c[i,j]*x[i,j],j=N minus {i}),i=N): 
conx:=seq( add(x[i,j], i=N minus {j})=1,j=N),  seq( add(x[i,j], j=N minus {i})=1,i=N):
conu:= seq(seq(u[i]-u[j]+n*x[i,j] <= n-1, i=N minus {1,j}), j=N minus {1}):

r := Optimization[LPSolve](z, {conx, conu}, #integervariables={seq(u[i],i=N minus {1})},
                           binaryvariables={ seq(seq(x[i,j],i=N minus {j}),j=N)} );

X:=eval(Matrix(n, symbol=x),{r[2][], seq(x[i,i]=0,i=1..n)}):
ord:=1: k:=1: to n do k:=max[index](X[k]); ord:=ord,k od:'order'=ord;

r := [22.6135858310497, [u[2] = -.999999998967403, u[3] = -2.99999997900383, u[4] = -8.88178419700125*10^(-16), u[5] = -2.00000000722817, u[6] = -3.99999998106902, x[1, 2] = 0, x[1, 3] = 0, x[1, 4] = 0, x[1, 5] = 0, x[1, 6] = 1, x[2, 1] = 0, x[2, 3] = 0, x[2, 4] = 1, x[2, 5] = 0, x[2, 6] = 0, x[3, 1] = 0, x[3, 2] = 0, x[3, 4] = 0, x[3, 5] = 1, x[3, 6] = 0, x[4, 1] = 1, x[4, 2] = 0, x[4, 3] = 0, x[4, 5] = 0, x[4, 6] = 0, x[5, 1] = 0, x[5, 2] = 1, x[5, 3] = 0, x[5, 4] = 0, x[5, 6] = 0, x[6, 1] = 0, x[6, 2] = 0, x[6, 3] = 1, x[6, 4] = 0, x[6, 5] = 0]]

order = (1, 6, 3, 5, 2, 4, 1)

If you need to read the same data several times you could do this:

1. Import once using ImportMatrix into A
2.    save A, "hugeA.m";
3.  Whenever you need A, use
     read "hugeA.m";
     It will be executed much faster (being in internal format .m).

The result is wrong, the implicitplot algorithm does not work well here.




Probably some "clever" options could be added, but always a better idea is to try a parametric representation of the curve:

restart;
f:=(1/2)*x*tan(log((1/5)*(2*(x^2+y^2)))) = y:
g:=simplify(eval(f, [x=r*cos(t), y=r*sin(t)])) assuming r>0:
u:=solve(g,t);

plot( [[r*cos(u),r*sin(u), r=0..5],[-r*cos(u),-r*sin(u), r=0..5]]);

The idea is similar to Kitonum's but has some advantages.
Note that some details will be omitted.

As it was observed,

f:=a^2+b^2+c^2+a*b+a*c+b*c-(a+b-c)*sqrt(2*a*b+a*c+b*c)-(a+c-b)*sqrt(a*b+2*a*c+b*c)-(b+c-a)*sqrt(a*b+a*c+2*b*c);
is positively homogeneous and symmetric, so we may suppose:
-1 <= a, b, c <= 1.

It is easy to see that c=0 implies a=b=0.
It will be enough to consoder c=1 and c=-1.

1.     c=1

f1:= eval(f, c=1);
plots:-implicitplot(f1, a=-1..1,b=-1..1, grid=[1000,1000]);

So, we have two branches. But they are symmetric because f(a,b,1)=f(b,a,1).
So we may consider only a>0, b<0.
Actually a will be in a1..a2  where a1, a2 could be computed [here at least of one sqrt must be 0].
So, the corresponding b will be h(a) = RootOf(f1, b=-1/5 .. -1/3).
[Note that h(a) can be expressed as a RootOf of a polynomial]

2.    c=-1
There are no solutions.

Conclusion:
The roots are:  a=k*u, b=k*h(u), c=k,  for all u in [a1,a2] and k>=0
and all permutations of these.

Edit.
 

a1:=-1+sqrt(1+RootOf(_Z^8-8*_Z^7-22*_Z^6+228*_Z^5+17*_Z^4-624*_Z^3-352*_Z^2+256*_Z+256, index = 2));

a2:=1;

h:= a -> 
RootOf(
a^16-8*a^15*b-4*a^14*b^2-48*a^13*b^3+154*a^12*b^4+40*a^11*b^5-304*a^10*b^6-56*a^9*b^7+451*a^8*b^8-56*a^7*b^9-304*a^6*b^10+40*a^5*b^11+154*a^4*b^12-48*a^3*b^13-4*a^2*b^14-8*a*b^15+b^16-8*a^15-8*a^14*b-192*a^13*b^2-56*a^12*b^3+96*a^11*b^4+456*a^10*b^5+296*a^9*b^6-560*a^8*b^7-560*a^7*b^8+296*a^6*b^9+456*a^5*b^10+96*a^4*b^11-56*a^3*b^12-192*a^2*b^13-8*a*b^14-8*b^15-4*a^14-192*a^13*b-480*a^12*b^2-1104*a^11*b^3+1564*a^10*b^4+3040*a^9*b^5-404*a^8*b^6-4048*a^7*b^7-404*a^6*b^8+3040*a^5*b^9+1564*a^4*b^10-1104*a^3*b^11-480*a^2*b^12-192*a*b^13-4*b^14-48*a^13-56*a^12*b-1104*a^11*b^2+720*a^10*b^3+4616*a^9*b^4+3464*a^8*b^5-4864*a^7*b^6-4864*a^6*b^7+3464*a^5*b^8+4616*a^4*b^9+720*a^3*b^10-1104*a^2*b^11-56*a*b^12-48*b^13+154*a^12+96*a^11*b+1564*a^10*b^2+4616*a^9*b^3+6008*a^8*b^4-1792*a^7*b^5-7416*a^6*b^6-1792*a^5*b^7+6008*a^4*b^8+4616*a^3*b^9+1564*a^2*b^10+96*a*b^11+154*b^12+40*a^11+456*a^10*b+3040*a^9*b^2+3464*a^8*b^3-1792*a^7*b^4-3896*a^6*b^5-3896*a^5*b^6-1792*a^4*b^7+3464*a^3*b^8+3040*a^2*b^9+456*a*b^10+40*b^11-304*a^10+296*a^9*b-404*a^8*b^2-4864*a^7*b^3-7416*a^6*b^4-3896*a^5*b^5-7416*a^4*b^6-4864*a^3*b^7-404*a^2*b^8+296*a*b^9-304*b^10-56*a^9-560*a^8*b-4048*a^7*b^2-4864*a^6*b^3-1792*a^5*b^4-1792*a^4*b^5-4864*a^3*b^6-4048*a^2*b^7-560*a*b^8-56*b^9+451*a^8-560*a^7*b-404*a^6*b^2+3464*a^5*b^3+6008*a^4*b^4+3464*a^3*b^5-404*a^2*b^6-560*a*b^7+451*b^8-56*a^7+296*a^6*b+3040*a^5*b^2+4616*a^4*b^3+4616*a^3*b^4+3040*a^2*b^5+296*a*b^6-56*b^7-304*a^6+456*a^5*b+1564*a^4*b^2+720*a^3*b^3+1564*a^2*b^4+456*a*b^5-304*b^6+40*a^5+96*a^4*b-1104*a^3*b^2-1104*a^2*b^3+96*a*b^4+40*b^5+154*a^4-56*a^3*b-480*a^2*b^2-56*a*b^3+154*b^4-48*a^3-192*a^2*b-192*a*b^2-48*b^3-4*a^2-8*a*b-4*b^2-8*a-8*b+1,
b, -3/10 .. -27/100)   :

To plot the solution use:

plot3d([ [k*u, k*h(u), k], [k*u, k, k*h(u)], [k*h(u), k*u,  k],
[k*h(u), k, k*u], [k, k*h(u), k*u], [k, k*u, k*h(u)] ],  u=a1..a2, k=0..10);