Find a closed surface for which $\iint_S \textbf{F} \cdot d\vec{S}$ is negative. $\text{div}\textbf{F} = x^2+y^2+1$. Find a closed surface for which $$\iint_S \textbf{F} \cdot d\vec{S}$$ is negative or otherwise state why it's not possible.

$\textbf{My attempt}$:
If $S$ and $D$ are a domain such that $\textbf{F}$ is continuous and differentiable on $S \cup D$, then by the divergence theorem:
$$\iint_S \textbf{F} \cdot d\vec{S} = \iiint_D  x^2+y^2+1 \ dV $$
The integral on the right is always positive and therefore there exists no such surface $S$ for which the flux integral becomes negative. 
Can anyone confirm that my approach is correct? Also, what if $S \cup D$ has singularities? Can we use the extended form of the divergence theorem:
$$\iint_S \textbf{F} \cdot d\vec{S} = \iint_{S'} \textbf{F} \cdot d\vec{S} +\iiint_{D'}  x^2+y^2+1 \ dV $$
and conclude that there may exist such a surface for which $$\iint_S \textbf{F} \cdot d\vec{S} < 0$$
$S'$ encloses the singularity.
 A: As you noted, your approach to the solution of the problem is entirely correct as long as the divergence theorem holds for the domain $D$ you consider, so your question is basically equivalent to asking for what kind of domains the Gauss-Green (divergence) theorem is true: and I'll try to answer to your question by answering to the more general one about the structure of domains for which the Gauss-Green (divergence) theorem holds in three steps.

*

*The divergence theorem is true for all domains $D$ whose perimeter i.e. the measure $|\partial D|=|S|$ of their boundary is finite: these sets are called Caccioppoli sets or sets of finite perimeter, and for them the general divergence formula
$$
\int\limits_D\!\operatorname{div} \mathbf{F}\, \mathrm{d}V 
=\int\limits_{S} \mathbf{F}\,\cdot \mathrm{d}\vec{S} =\int\limits_{S} \mathbf{F}\,\cdot\nu_D\, \mathrm{d}|\nabla\chi_D|\quad\forall\mathbf{F}\in [C_c^1(\mathbb{R}^n)]^n\tag{1}\label{1}
$$
where $\chi_D$ is the characteristic (indicator) function of the domain $D$ and $|\nabla\chi_D|$ is the (finite) total variation measure associate to its (generalized) gradient, holds true. You can find a sketch of the proof and detailed references in this Q&A. Formula \eqref{1} holds irrespective of the kind of singularity $S=\partial D$ may have: there can be cusps, edges, or wedges alike. Thus, for any positive divergence vector field $\mathbf{F}$, the flux across the boundary of any set of finite perimeter is necessarily positive.


*However, there are also classes of sets which do not necessarily have a finite perimeter for which the divergence theorem holds. As pointed out by Ben McKay in his answer to this question, Friedrich Sauvigny  proves the general Stokes theorem for sets whose boundary has a singular part with finite capacity ([1], §1.4 pp. 30-39) and, as a corollary, a generalized Gauss-Green (divergence) theorem is demonstrated ([1], §1.5, pp. 39-49): and the class of finite capacity sets considered by Sauvigny strictly includes the class of Caccioppoli sets. Thus again, for any positive divergence vector field $\mathbf{F}$, the flux across any closed boundary whose singular part has a finite capacity (in the sense of [1]) is necessarily positive.


*Finally, it should be noted that the accepted answer of Paul Siegel to the question above, goes further by conjecturing that a form of Stokes theorem in the plane may hold for every proper subset: therefore, conjecturing that his "bold" assertion is true also for higher dimensional regions, it may be difficult if not impossible to find a set for which the flux across it of a positive divergence vector field $\mathbf{F}$ is negative, even if an infinite perimeter/capacity (in the sense of [1]) is allowed.
Final conclusions.
As you have noted in the premise to your question, if the Gauss-Green-Divergence theorem holds for a given domain $D$ then any positive divergence vector field $\mathbf{F}$ would generate a positive flux across its boundary. Summing up, if $\operatorname{div}\mathbf{F}=x^2+y^2+1$ then

*

*It is not possible to find a domain $D$ with a finite perimeter for which
$$
\int\limits_{S} \mathbf{F}\,\cdot \mathrm{d}\vec{S} <0\label{2}\tag{2}
$$
since for those domains (Caccioppoli sets) the Gauss-Green (divergence) formula \eqref{1} holds

*In the same way, it is not possible to find a domain $D$ whose singular part of the boundary has a finite capacity in the sense of [1] for which \eqref{2} holds, since even for those domains identity \eqref{1} holds.

*And even if you consider sets of infinite perimeter whose singular part  of their boundary has an infinite capacity, it may be nevertheless be impossible for \eqref{2} to hold, since it could be possible that, in a properly generalized sense, \eqref{1} may still hold.

Reference
[1] Sauvigny, Friedrich, Partial differential equations 1. Foundations and integral representations. With consideration of lectures by E. Heinz, 2nd revised and enlarged ed. (English), Universitext, Berlin: Springer Verlag, ISBN 978-1-4471-2980-6/pbk; 978-1-4471-2981-3/ebook, pp. xv+447 (2012), MR2907678, Zbl 1246.35001.
