Integration over arbitrary domains In mathematical physics, we sometimes encounter situations where we have integrals of the forms:
$$\text{(case 1):}\ \ \ \ \int\limits_{D} f(x,y,z) dx dy dz =k$$
 $$\text{(case 2):}\ \ \ \ \int\limits_{D} f(x,y,z) dx dy dz =\int\limits_{D} g(x,y,z) dx dy dz$$
where $D$ is an arbitrary domain (e.g. volume) of intergation over the $(x,y,z)$ variables and $k$ is some constant. 
Because of $D$'s arbitrariness (i.e. can be taken small enough so that the integrand is taken outside integral sign), we usually proceed by taking $f,g$ outside the integral and ending up with the corresponding results:
$$\text{(case 1):}\ \ \ \ f=k/V_{D}$$
 $$\text{(case 2):}\ \ \ \ f=g$$
where $V_{D}$ is the volume $\int_{D}dxdydz$.
Questions:
(a) Is there any formal/rigorous theorems in math (e.g. in analysis) that specialise in such situations or explain its conditions in general? What branch/subbranch of math would cover such problems?  
(b) If we are faced with a problem of case (2) above, for example, and we could find some kind of domains $D$ of a certain shape but arbitrary size (e.g. say $D$ as a sphere of arbitrary radius) for which both sides (integrals) are known to be equal (and hence we reach equality of integrands), will that be enough to conclude that equality of integrals will reduce to equality of integrands in general (or do we need to first show that it is also true for all other arbitrary shapes of $D$)? 
(c) Finally, could the same discussion/conclusions of case (2) above be extended to the case of integrals over closed surfaces (via Divergence's theorem for instance?): $\oint\limits_{A} \boldsymbol{F}\cdot \boldsymbol{dA}=\oint\limits_{A} \boldsymbol{G}\cdot \boldsymbol{dA}$?
(Note - I have asked this question on math SE, but got no answers.) 
 A: You may be interested in the book "Methods of Applied Mathematics" by Arbogast and Bona, available here as a pdf 
https://www.google.co.uk/url?q=https://www.ma.utexas.edu/users/arbogast/appMath08c.pdf&sa=U&ved=0ahUKEwjL5s2U-MPVAhWjLMAKHenGAAMQFggOMAA&usg=AFQjCNGW21X85j7nqyk6zuDOxRG5wdZ5Pg
Its quite a tough read, but of interest is proposition 1.39 on page 21, which can be restated for continuous functions that, if for all measurable $D \subset \Omega$, where $\Omega$ is the domain of our function (and a union of finite balls, so that we dont have any wierd measure zero bits), we have
$$
\int_D f dV = 0
$$
Then
$$
f(x) = 0 \forall x \in \Omega
$$
This does you for case 1 when the right hand side is zero. I'm not at all sure that the nonzero case is even true...
For case 2
$$
\int_D f dV = \int_D g dV
$$
$$
\int_D f - g dV = 0
$$
$$
f-g=0
$$
$$
f=g
$$
I hope that answers part a). For part b) the answer is no. If you have a series of concentric spheres for which the volume integrals are equal, then you can then conclude that the surface integrals over the spheres are equal, but the values may be distributed differently over those spherical surfaces. 
For part c), if those surfaces are not necessarily closed, then a version the above theorum will apply and you have that $F=G$ (I dont have a proof of this, just intuition). If the surfaces must be closed then the theorum will not apply since we do not have the result for all domains, but we can use the divergence theorum and then the above theorum to obtain that $\nabla \cdot F=\nabla \cdot G$
Edit
Here I correct and expand on my answer to part b). 
For a simply connected open set $D$ and a continuous function $f$ we have the multivariate mean value theorem for integrals, which states that there exists some $c \in D$ for which
$$
\int_D f(x) dV = f(c) |D|
$$
where $|D|$ is the volume (specifically, the measure) of the open set. 
Therefore, if we have an infinite sequence of simply connected open sets $D_0,D_1,...$ such that the sets all include some point $x_0$, and the sets 'tend towards' this point, specifically
$$
\lim_{n \rightarrow \infty} (\max_{x \in D_n} |x-x_0|) = 0
$$
then for each domain we can specify a point $c_n$ for which
$$
\frac{1} {|D_n|} \int_{D_n} f(x) dV = f(c_n) =0
$$
where the equality with zero is from the question. We have that $c_n \in D_n$ thus $c_n \rightarrow x_0$ as $n \rightarrow \infty$ and $f(x_0)=0$.
Thus, if you find a set of simply connected open domains of arbitrarily small size around every point, the above shows that $f=0$.
We can therefore state the following theorem:
If, for each point $x_0$ in an open domain $\Omega$, we have that, for a sequence of simply connected open subsets $D_n$, $n$ any natural number, 
$$
\int_{D_n} f(x) dV =0
$$
for some continuous function $f$, and
$$
\lim_{n \rightarrow \infty} (\max_{x \in D_n} |x-x_0|) = 0
$$
then
$$
f=0
$$
A: I am here posting a second answer, which will be a complete restructuring of the first, inclusion of discussion in the comments, and expansion of the scope of the answer. I will leave it to others to decide whether to keep or delete the other.
Firstly, it is important that "case 1" is only valid for $k=0$. "Case 2" is equivalent to it because
$$
\int_D f dV = \int_D g dV
$$
$$
\int_D (f - g) dV = 0
$$
then if "case 1" is true
$$
f-g=0
$$
$$
f=g
$$
therefore "case 1" implies "case 2".
Part c) of the question is also equivalent to "case 1", by the divergence theorem (as noted in the question), if $A$ is the (piecewise smooth) boundary of the bounded open set $D$ then
$$
\oint_A f \cdot \hat{n} dA = \int_D \nabla \cdot f dV
$$
thus if "case 1" is true and the and the integrals over all (or sufficiently many) surfaces are zero then 
$$
\nabla \cdot f =0
$$
However, if the surfaces are not necessarily closed then, whilst the above argument holds, a stronger case can be made and we will find that $f=0$. Note that this implies that $\nabla \cdot f =0$ and there is no contradiction.

Our first job, then, is to argue "case 1". This is a standard result in the field of functional analysis, see Methods of Applied Mathematics by Arbogast and Bona, available here as a pdf. Its quite a tough read, but of interest is proposition 1.39 on page 21, which can be restated for continuous functions.

Theorem 1 Let f be a continuous function whose domain is an open set $\Omega$. If for all measurable $D \subset \Omega$ we have
  $$
\int_D f dV = 0
$$
  then
  $$
f(x) = 0 \forall x \in \Omega
$$

Note that 'measurable' just means that a Lebesgue volume integral can be defined over such a set, and so this is a very strong requirement on the function. In line with part b) of the question, we can construct a slightly more complicated theorem, which requires us to show the 'integral=0' condition on a smaller collection of sets.

Theorem 2 If, for each point $x_0$ in an open domain $\Omega$, we have that, for a sequence of simply connected open subsets $D_n$, $n$ any natural number, 
  $$
\int_{D_n} f(x) dV =0
$$
  for some continuous function $f$, and
  $$
\lim_{n \rightarrow \infty} (\max_{x \in D_n} |x-x_0|) = 0
$$
  then
  $$
f(x) = 0 \forall x \in \Omega
$$
Proof For a simply connected open set $D$ and a continuous function $f$ we have the multivariate mean value theorem for integrals, which states that there exists some $c \in D$ for which
  $$
\int_D f(x) dV = f(c) |D|
$$
  where $|D|$ is the volume (specifically, the Lebesgue measure) of the open set.
  Therefore, if we have an infinite sequence of simply connected open sets $D_0,D_1,...$ such that the sets all include some point $x_0$, and the sets 'tend towards' this point, specifically
  $$
\lim_{n \rightarrow \infty} (\max_{x \in D_n} |x-x_0|) = 0
$$
  then for each domain we can specify a point $c_n$ for which
  $$
\frac{1} {|D_n|} \int_{D_n} f(x) dV = f(c_n) =0
$$
  where the equality with zero is because we have asserted that the integral is zero. We have that $c_n \in D_n$ thus $c_n \rightarrow x_0$ as $n \rightarrow \infty$ and $f(x_0)=0$.
  Thus, if you find a set of simply connected open domains of arbitrarily small size around every point, the above shows that $f=0$.

This covers question a) and also question b), it is fine to use open domains of a particular shape and arbitrary size, so long as they can be centred around any point and limited to zero size. We could also use different shapes at different sizes, so long as the conditions of the above are met.

At this point I will sidestep slightly from analysis, and address continuum mechanics, and the approximations of mathematical physics. In physical systems there is typically a split between the statistical macroscopic properties and the exact models that must be used to capture the small scale motion. The small scale motion may be the vibration of molecules in a body of water, or the bouncing of stones in an avalanche. We typically see that the statistical macroscopic variables accurately capture the large-scale behaviour, and the details of the microscopic are unimportant.
In such a case we typically have three scales of space-time. A large scale over which the dynamics of interest occurs so that the statistical variables are rapidly varying but approximately integral=0. A mid scale over which the statistical variables are slowly varying and approximately integral=0. A small scale at which the integral is no longer valid because the small scale behaviour dominates. Note that the statistical variables will be averages over a region of space-time of the same scale as the smallest regions at which the integral equation is valid.
What is typically then done is to define a macroscopic approximation to the physical system, where we claim that, since we have
$$
\int_{D_n} f(x) dV \approx 0
$$
for a sequence of domains satisfying
$$
\lim_{n \rightarrow \infty} (\max_{x \in D_n} |x-x_0|) = L
$$
where $L$ is the mid space-time scale (much smaller than the length scale of the dynamics), then we can approximate by
$$
f(x) = 0 \forall x \in \Omega
$$
Notice how similar this is to Theorem 2. It is not something that can be 'proved' mathematically (or at least, not that I have ever seen...), rather an approximation that often works.
Finally, it is worth pointing out here that none of this discussion relied on $\Omega$ being a subset of three dimensional Euclidean space. The open domain of $f$ could be, say, a surface of a fluid, or any other 'boundary' of the physical domain which will be, in itself, a physical domain. This is how we can think of boundary conditions, down to a very small area (or length, etc) of a boundary domain, the equation integral=0 holds to a good approximation. Of course, to get any of the micro-scale behaviour the integral will have to have some thickness, but the length scale across the boundary will be much smaller than the length scale along the boundary (this is what makes it a boundary), so the domains of integration will be very thin disks. We can use these in Theorem 2 (ignoring the thin thickness and the fact that the limit can only be evaluated down to some small length scale) which yields
$$
f=0
$$
on the 'boundary'.
A: Case 2 is evident, if it's satisfied for EVERY possible domain $D$.
As for case 1, it looks logical too, but if you're looking for a rigurous way, maybe you can do it using the Dirac Delta Function to re-write your domain. The properties of the Diract Delta function yield your result. 
