Here is a simple volume.

MyVolume=RegionIntersection[
Ball[{0,0,0},2], 
Cuboid[{0,0,0},{2,2,2}]
];
MyPic=BoundaryDiscretizeRegion[MyVolume,
Axes->True,
Boxed->True,
PlotTheme->"Lines",
BoxRatios->Automatic,
AxesLabel->{"x","y","z"}
]

It is easy to show the outward normals on it.  Note: they are all the same length and all outwards.

BoundaryDiscretizeRegion[MyVolume,
PlotTheme->"FaceNormals",
Boxed->True]

Integrating over the volume is easy.

g[{x_,y_,z_}]:= {-(x^2)+z,ArcTan[y+z+x]+ Sin[(E^(-x+(z^3)))],z+ Sin[z+ x y] }
Divg[{x_,y_,z_}]=Div[g[{x,y,z}],{x,y,z}]

NIntegrate[Divg[{x,y,z}],
{x,0,2},
{y, 0, (Sqrt[4-(x^2)])},
{z,0, (Sqrt[4-(x^2)-(y^2)])}]

NIntegrate[Divg[{x,y,z}], {x,y,z}∈MyVolume]

The flux integral comes in four pieces: three flat bits and a bit of a sphere.  
Each of the flat bits is a piece of a circle:  polar coordinates in the appropriate pair of variables for those;  the normals are {-1,0,0},{0,-1,0},{0,0,-1}.
Use spherical coordinates with ρ=1 and angle ranges 0<ϕ<π/2 and 0<θ<π/2 for the bit of the sphere; the normal is {ρ Cos[θ]Sin[ϕ],ρ Sin[θ]Sin[ϕ], ρ Cos[ϕ]}.

ρ=2;
Bits={NIntegrate[
g[{ρ Cos[θ]Sin[ϕ],ρ Sin[θ]Sin[ϕ],ρ Cos[ϕ]}].{ρ Cos[θ]Sin[ϕ],ρ Sin[θ]Sin[ϕ],ρ Cos[ϕ]} ρ Sin[ϕ],
{θ,0, π/2},
{ϕ,0, π/2}
],
NIntegrate[ g[{x,y,0}].{0,0,-1},
{x,0,2},
{y,0,(Sqrt[4-(x^2)])}], 
NIntegrate[ g[{x,0,z}].{0,-1,0},
{x,0,2},
{z,0,(Sqrt[4-(x^2)])}], 
NIntegrate[ g[{0,y,z}].{-1,0,0},
{y,0,2},
{z,0,(Sqrt[4-(y^2)])}]
} 
Total[Bits]

MyMeshPoints=RandomSample[
MeshCoordinates[BoundaryDiscretizeRegion[MyVolume]],
200];
Show[
MyPic,
VectorPlot3D[ 
g[{x,y,z}],
{x,0,2},
{y,0,2},
{z,0,2},
VectorStyle->Pink,
VectorPoints->MyMeshPoints,
BoxRatios->Automatic
]
]

Here is a more interesting volume.

MyVolume=ImplicitRegion[(x^2+(9/4)y^2+z^2-1)^3-x^2z^3-(9/80)y^2z^3≤0,{x,y,z}];
BoundaryDiscretizeRegion[MyVolume,
Axes->True,
Boxed->True,
BoxRatios->Automatic,
AxesLabel->{"x","y","z"}
]

It is easy to show the outward normals on it.  Note: they are all the same length and all outwards.

BoundaryDiscretizeRegion[MyVolume,
PlotTheme->"FaceNormals"]

Integrating over the volume is still easy. The flux integral would be MUCH harder.

g[{x_,y_,z_}]:= {(x^2),y,z};
Divg[{x_,y_,z_}]=Div[g[{x,y,z}],{x,y,z}]
NIntegrate[Divg[{x,y,z}], {x,y,z}∈MyVolume]//Chop

Break

DiscretizeRegion[ Polygon[{{0,0},{2,0},{2,1},{0,1}}],
Method->"DiscretizeGraphics"]
DiscretizeRegion[ Polygon[{{0,0},{2,0},{2,1},{0,1}}],
Method->"RegionPlot"]
DiscretizeRegion[ Polygon[{{0,0},{2,0},{2,1},{0,1}}],
Method->"Continuation"]
DiscretizeRegion[ Polygon[{{0,0},{2,0},{2,1},{0,1}}],
Method->"Boolean"]