> restart:
> with(linalg):
> phi:=1/sqrt(x^2+y^2+(z-1)^2)-1/sqrt(x^2+y^2+(z+1)^2);
                      1                               1              
        ----------------------------- - -----------------------------
                                (1/2)                           (1/2)
        / 2    2    2          \        / 2    2    2          \     
        \x  + y  + z  - 2 z + 1/        \x  + y  + z  + 2 z + 1/     
> u:=grad(-phi,[x,y,z]);n:=vector(3,[0,0,-1]);


> un:=evalm(transpose(n)&*u);
                   2 z - 2                           2 z + 2            
     - ------------------------------- + -------------------------------
                                 (3/2)                             (3/2)
         / 2    2    2          \          / 2    2    2          \     
       2 \x  + y  + z  - 2 z + 1/        2 \x  + y  + z  + 2 z + 1/     
> f:=subs(z=0,un);
                                     2         
                             ------------------
                                          (3/2)
                             / 2    2    \     
                             \x  + y  + 1/     
> assume(y>0,x>0);
> sol:=int(int(f,x=0..(1-y/2)),y=0..1);
                                  /1  (1/2)\   1   
                         -2 arctan|- 2     | + - Pi
                                  \4       /   2   
> evalf(sol);
                                0.8911225082
> tau:=vector(2,[(1-xi)*(1-eta)*(i-1)/N*(1-(j-1)/2/N)+(xi)*(1-eta)*(i)/N*(1-(j-1)/2/N)+(xi)*(eta)*(i)/N*(1-(j)/2/N)+(1-xi)*(eta)*(i-1)/N*(1-(j)/2/N),(1-xi)*(1-eta)*(j-1)/N+(xi)*(1-eta)*(j-1)/N+(xi)*(eta)*(j)/N+(1-xi)*(eta)*(j)/N]);

> taup:=map(simplify,jacobian(tau,[xi,eta]));tau:=map(simplify,tau);


> d:=det(taup);
                              2 N - j + 1 - eta
                              -----------------
                                       3       
                                    2 N        
> N:=3;
                                      3
> om:=vector(3,[5./18,4./9,5./18]);gp:=vector(3,[1./2-sqrt(15.)/10,1/2,1./2+sqrt(15.)/10]);


> sumand:=om[k]*om[l]*subs(x=subs(xi=gp[k],eta=gp[l],tau[1]),y=subs(xi=gp[k],eta=gp[l],tau[2]),f*subs(eta=gp[l],d));
    /              /7    1      1       \\///                            
    |2 om[k] om[l] |-- - -- j - -- gp[l]|| ||/7      1      7    7       
    \              \54   54     54      // |||-- i + -- j - -- + -- gp[k]
                                           \\\18     18     18   18      

         1          1                1            1        1         \  
       + -- gp[l] - -- gp[k] gp[l] - -- gp[k] j - -- i j - -- gp[l] i|^2
         18         18               18           18       18        /  

                            2    \      \
         /1     1   1      \     |      |
       + |- j - - + - gp[l]|  + 1|^(3/2)|
         \3     3   3      /     /      /
> evalf(sum(sum(sum(sum(sumand,k=1..3),l=1..3),i=1..N),j=1..N)-sol);
                                         -8
                                  -4.6 10  
>