Integrating a function over a surface integral

Question

From a first principles bandstructure calculation I get an energy scalar field in three dimensions $E(x,y,z)$. It's now easy to plot a constant energy (contour)-surface for dedicated values $\epsilon=E(x,y,z)$ by using ListContourPlot3D. Lets assume the following example scalar field (e.g. the case of a graviatation potential field):

data = Table[1/Sqrt[x^2 + y^2 + z^2], {x, -1, 1, 0.1}, {y, -1, 1, 0.1}, {z, -1, 1,0.1}]
        /. ComplexInfinity -> 10^16 // Quiet;

I exchanged the infinite value at $x=y=z=0$ with a finite one to prevent any error messages coming from the visualisation routines. Now the iso-surfaces can be easily plotted for e.g. $\epsilon=2.0$

ListContourPlot3D[data, Contours -> {2.0}, DataRange -> {{-1, 1}, {-1, 1}, {-1, 1}}]

Coming back to the bandstructure calculation, the energy surfaces are much more complicated and cannot be described analytically for the general case as can be seen from the following plot of an energy surface for Silicon:

An important electronic property, the DOS (Density of States) is calculated through the integral $$\oint \limits_{E(x,y,z)=\epsilon}\frac{dS}{\lvert\nabla E(x,y,z)\rvert}.$$

and additional integrals of interest have the general form

$$\oint \limits_{E(x,y,z)=\epsilon}\frac{dS}{\lvert\nabla(\Delta E(x,y,z))\rvert}A(x,y,z)$$

with $A(x,y,z)$ as some additional scalar function (e.g. for effective mass tensor, $A(x,y,z)$ is the tensor product of the gradients of the energy).

Now for the question. How to most efficiently implement such integrals on numerically specified data (see above) over iso-surfaces in Mathematica?

EDIT: Trying out the proposal by Mastrok I run into convergence problems for e.g. $\epsilon$ = 2.0 of the integration (most probably due to the singularity in the coordinate center and the additional derivatives on the interpolation function when using Gauss' Theorem):

e = ListInterpolation[data, {{-1, 1}, {-1, 1}, {-1, 1}}];
dive = Simplify[Div[#/Norm[#]^2& @ Grad[e[x, y, z], {x, y, z}], {x, y, z}] /. 
  Derivative[1][Abs][x_] :> x/Abs[x] /. {x->#1, y->#2, z->#3} &];
NIntegrate[Piecewise[{{dive[x, y, z], e[x, y, z] < 2.0}}, 0] ,
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1}
]

(*Output*)
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.

Answer 1

First note that,

$$ \oint_{E(\vec{x})=\epsilon}\frac{dS}{|\nabla E(\vec{x})|} =\oint_{E(\vec{x})=\epsilon}\frac{\nabla E(\vec{x})\cdot d\vec{S}}{|\nabla E(\vec{x})|^2} = \int_{E(\vec{x})\leq\epsilon} \nabla\cdot\left( \frac{\nabla E(\vec{x})}{|\nabla E(\vec{x})|^2}\right) dV$$

By now, the volume integral can be evaluated by

Integrate[f Boole[E[x,y,z]<epsilon],{x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax}]

where $f$ is the integrant and

 Boole[E[x,y,z]<epsilon] 

forces the integrant to be 0 outside the region. Therefore, you need to have large enough {x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax} to cover your whole domain.

Answer 2

By using the output of ListContourPlot3D and the new Mathematica 10.0 feature DiscretizeGraphics, one can nicely generate a meshed contour region which is suitable for NIntegrate. We can show this for the above example for energy contours from 1.0 to 2.0 with a step width of 0.1:

Monitor[Table[
  e = ListInterpolation[data, {{-1, 1}, {-1, 1}, {-1, 1}}];
  f = Simplify[1/Norm[Grad[e[x, y, z], {x, y, z}]]];
  plot = ListContourPlot3D[data, Contours -> {ϵ}, 
   DataRange -> {{-1, 1}, {-1, 1}, {-1, 1}}];
  R = DiscretizeGraphics[Normal[plot /. (Lighting -> _) :> Lighting -> Automatic]];
  NIntegrate[f, {x, y, z} ∈ R], ϵ, 1.0, 2.0, 0.1}], ϵ]
, {ϵ, 1.0, 2.0, 0.1}], ϵ]
(*Output*)
{12.5221, 8.55496, 6.03428, 4.38027, 3.25017, 2.46631, 1.90372, 1.49193, 1.18722, 0.952876, 0.777326}

The result is very close to the analytical solution. The differences are coming from the low sampling rate of the original data generation:

anasol = Table[
 Integrate[
  Simplify[1/Norm[Grad[1/Sqrt[x^2 + y^2 + z^2], {x, y, z}]]], 
   {x, y, z} ∈ 
    ImplicitRegion[x^2 + y^2 + z^2 == (1/ϵ)^2, {x, y, z}]],
  {ϵ, 1.0, 2.0, 0.1}]
(*Output*)
{12.5664, 8.583, 6.06017, 4.39984, 3.27113, 2.48225, 1.91748, 1.50458, 1.19707, 0.964263, 0.785398}

Answer 3

I've been try to perform calculations of exactly this type, and trying the different methods suggested. The energy surface I need is just as convoluted as the silicon Fermi surface shown in the OP and highly computationally intensive, and the scalar function I am trying to calculate the integral of over the surface (the function called $A(x,y,z)$ in the OP) is also very computationally intensive. Finding a mesh for the surface, and the gradient were too much, and there isn't enough control provided in the available routines. So in my case none of the methods above have been able to work with anything like a reasonable convergence in an acceptable computation time. The best method for speed that I've found is kind of a brute-force method, and it produces the entire spectrum in a smaller timescale than the surface integral method for a single energy (both in my case, and for the OP minimal code data). As suggested in earlier replies, one can turn the surface integral into a volume integral. Unlike the volume integral version suggested by Mastrok, which still requires finding the surface and gradient (which was far too computationally demanding for me), a much simpler volume integral is produced by $$ \int_{E(x,y,z)=\epsilon}dS \frac{A(x,y,z)}{|\nabla E(x,y,z)|}= \int A(x,y,z)\delta(E(x,y,z)-\epsilon) dV $$ where $\delta(E(x,y,z)-\epsilon)$ is a Dirac delta function. The key here is to allow the delta function to have a finite width in energy, and remember that it has unit area so that it has a height equal to the reciprocal of the width. Then you can just evaluate A and E on a grid, and then bin and sum the A values according to the E's (treating x,y,z as dummy parameters).

r = 0.1; (* resolution in x,y,z *)
data = 
 Table[1/Sqrt[x^2 + y^2 + z^2], {x, -1, 1, r}, {y, -1, 1, r}, {z, -1, 
     1, r}] /. ComplexInfinity -> 10^16 // Quiet; 
dataA = 
 Table[1, {x, -1, 1, r}, {y, -1, 1, r}, {z, -1, 1, r}];
ew = 0.1; (* energy width of the delta function *)
surfint = HistogramList[
   WeightedData[Flatten[data], Flatten[dataA]], {1 - ew/2, 2 + ew/2, 
    ew}];
surfint = surfint[[2]] r^3/ew; (* take care of the resolution of the grid and the delta function width *)
ListPlot[{surfint, anasol}, DataRange -> {1.0, 2.0}]

The above code produces the following (where anasol, shown in yellow, is the exact result from another earlier answer). Note that in this example dataA is just a grid of 1's so that the result could be compared with the exact solution, but it could be any function of x,y,z. The result is a little more noisy than other answers, but far faster, and so the noise can be improved by reducing the resolution and playing with bin widths.