VC.10 3D Surface Measurements
 BASICS

\(Off[General::spell]; 
Off[General::spell1]; 

    \(Off[ParametricPlot3D::ppcom]; \)\)Off[ParametricPlot::ppcom]; Off[Plot::plnr]; SetOptions[Limit, Analytic -> True];

CMView = {2.7, 1.6, 1.2};

Vector Drawer

BeginPackage["Vector3D`", "Graphics`Colors`"];

\(Axes3D::"usage" = 
      "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8]."; 

    \(Perpend::"usage" = 
      "Perpend[a] returns a unit vector which is perpendicular to the vector a."; \)\)Arrow::"usage" = 
      "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. \
Arrow[a, Tail -> tail] gives a vector from tail to a."; Vector::"usage" = 
      "Vector[a, b] gives a 2D or 3D vector from a to b."; ArrowHead::"usage" = 
      "ArrowHead[a, vec] produces an arrowhead with its tip placed at the \
point a, pointing in the direction of the vector vec."; Tail::"usage" = 
      "Tail -> point puts the tail of the vector at point."; Norm::"usage" = 
      "Norm[vec] returns the Euclidean length of vec."; Normalize::"usage" = 
      "Normalize[vec] returns a unit vector in the direction of vec, \
provided vec is nonzero."; Aperture::"usage" = 
      "Aperture is the ratio of the radius of the base to the length of the \
head of a vector."; TipSize::"usage" = 
      "TipSize specifies an absolute size for the inner tip length of the \
head of a vector."; TipRatio::"usage" = 
      "TipRatio is the ratio of the inner length of the arrowhead to the \
length of the shaft of a vector."; HeadRatio::"usage" = 
      "HeadRatio is the ratio of the outer length of the arrowhead to the \
length of the shaft of a vector."; HeadSize::"usage" = 
      "HeadSize specifies an absolute size for the head of a vector."; 
    ScaleFactor::"usage" = 
      "ScaleFactor specifies the amount to scale a vector in length and may \
be either a number or a pure function."; ZeroVectorPointSize::"usage" = 
      "ZeroVectorPointSize is the size of the point used to represent a \
zero vector."; Options[Vector] = 
      Options[Arrow] = 
        Options[ArrowHead] = {HeadRatio -> 0.18, 
            TipRatio -> 0.14, Aperture -> 0.3, 
            ScaleFactor -> 1, HeadType -> Polygon, 
            ShaftQ -> True, EdgesQ -> True, 
            VectorColor -> RGBColor[0, 0, 1], ShaftWidth -> 0.005, 
            ZeroVectorPointSize -> 0.01}; Begin["`Private`"]; SetOptions[ParametricPlot, AspectRatio -> Automatic]; SetOptions[Plot, AspectRatio -> Automatic]; SetOptions[Graphics, AspectRatio -> Automatic]; Axes3D[u_, v_] := 
      Graphics3D[{{Blue, Line[{{-(u/3), 0, 0}, {u, 0, 0}}]}, 
          Text["x", {u + v, 0, 0}], {Blue, 
            Line[{{0, -(u/3), 0}, {0, u, 0}}]}, 
          Text["y", {0, u + v, 0}], {Blue, 
            Line[{{0, 0, -(u/3)}, {0, 0, u}}]}, 
          Text["z", {0, 0, u + v}]}]; Axes3D[u_] := Axes3D[u, (u/8)]; Norm[a_?VectorQ] := N[(Sqrt[a . a])]; Normalize[a_?VectorQ] := If[a . a == 0, a, (a/Norm)[a]]; Perpend[a_?VectorQ] := 
      Normalize[{a[[2]], 
            -a[[1]]}] /; 
        Length[a] == 2Perpend[a_?VectorQ] := 
      Normalize[
          If[a[[1]] == 0, {1, 0, 0}, {0, 
              a[[3]], 
              -a[[2]]}]] /; 
        Length[a] == 3Base3D = 
      Table[{Cos[((2 k π)/8.) ], Sin[((2 k π)/8.) ]}, {k, 0, 
          8}]; Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = 
            {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, 
                  TipRatio, Aperture, ScaleFactor, VectorColor, 
                  ZeroVectorPointSize} /. {opts} /. Options[Vector]; 
          If[from == to, 
            Return[Graphics[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]; tr = 0.75 hr]; 
          If[NumberQ[ts], tr = (ts/Norm)[diff]]; 
          Graphics[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              With[{trans = ap hr Norm[diff] Perpend[diff]}, {vc, 
                  Thickness[sw], 
                  ht[{tip, tip - hr diff + trans, tip - tr diff, 
                      tip - hr diff - trans, tip}]}]}]] /; 
        Length[from] == Length[to] == 2Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, eq, sw, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = 
            {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, 
                  ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts}
                   /. Options[Vector]; 
          If[from == to, 
            Return[Graphics3D[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]]; tr = (hr/2); 
          Graphics3D[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              If[eq, {}, EdgeForm[]], {SurfaceColor[vc], 
                With[{temp = Perpend[diff]}, 
                  (Polygon[Append[#1, tip]]&)/@
                    Partition[
                      (tip - hr diff + #1&)/@
                        ((ap hr Norm[diff] Base3D) . {temp, 
                              Normalize[temp\[Cross]diff]}), 2, 1]]}}]] /; 
        Length[from] == Length[to] == 3Arrow[a_, Tail -> b_, opts___] := Vector[b, a + b, opts]Arrow[a_, opts___] := Vector[Table[0, {Length[a]}], a, opts]\(ArrowHead[a_, b_, opts___] := 
      Vector[a - b, a, ShaftQ -> False, opts]; 

End[]; 

    \(EndPackage[]; \)\)

Graphics`Colors`GosiaGreen = RGBColor[0,  0.392187,  0];

Spherical Setup

\(Clear[polar]; 

    polar[{x_, y_, z_}] := {\@\(x\^2 + y\^2 + z\^2\), ArcTan[y\/x], 
        ArcCos[z\/\@\(x\^2 + y\^2 + z\^2\)]}; 
Clear[X, length]; 

    length[X_] := \@\(X . X\); 
\(rr = 3; \)\)R = 5; scale = 0.8; red = RGBColor[1, 0, 0]; blue = RGBColor[0, 0, 1]; color1 = RGBColor[0.8, 0.2, 0.8]; color2 = RGBColor[0.2, 0.8, 0.2]; black = RGBColor[0, 0, 0]; \(Clear[axes]; 

    axes := Graphics3D[{Line[{X, {1.5\ rr, 0, 0}}], 
          Line[{X, {0, 1.5\ rr, 0}}], Line[{X, {0, 0, 2\ rr}}]}]; 

    Clear[xyplane]; 

    xyplane := 
      Graphics3D[
        Polygon[{{\(-2\)\ rr, \(-2\)\ rr, 0}, {\(-2\)\ rr, 2\ rr, 0}, {2\ rr, 
              2\ rr, 0}, {2\ rr, \(-2\)\ rr, 0}}]]; 
\(X = {0, 0, 0}; \)\)origin = Graphics3D[{black, Point[X]}]; Z = {2, 3, 5}; \(Clear[point]; 

    point := Graphics3D[{red, PointSize[0.03], Point[scale\ Z]}]; 

    \(Z0 = {0, 0, length[0.8\ Z]}; \)\)\(Clear[zpoint]; 
zpoint = Graphics3D[{black, Point[Z0]}]; 

    \(XY = scale\ Z\ {1, 1, 0}; \)\)\(Clear[xypoint]; 
xypoint := Graphics3D[{black, Point[XY]}]; 

    \(X0 = {length[scale\ Z\ {1, 1, 0}], 0, 0}; \)\)xpoint = Graphics3D[{black, Point[X0]}]; Y0 = {0, length[scale Z {1, 1, 0}], 0}; ypoint = Graphics3D[{black, Point[Y0]}]; \(Clear[ray]; 
ray := Graphics3D[{blue, Line[{X, scale\ Z}]}]; 

    Clear[ray1]; 
ray1 := Graphics3D[{color1, Line[{scale\ Z, Z}]}]; 

    Clear[xyray]; 
xyray := Graphics3D[{red, Line[{X, Z\ {1, 1, 0}}]}]; 

    \(rays := {ray, xyray}; \)\)xyaux = 
      (length[0.8 Z] (0.8 Z {1, 1, 0}))/
        length[0.8 Z {1, 1, 0}]; \(Clear[zarc]; 

    zarc = Graphics3D[{blue, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 0, 
                \(polar[0.8\ Z]\)[[3
                  ]],  .1}]]}]; 
Clear[zarc1]; 

    zarc1 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(polar[0.8\ Z]\)[[3]], 
                1,  .1}]]}]; 
Clear[zarc2]; 

    zarc2 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(-0.8\), 0,  .1}]]}]; 
Clear[xyarc]; 

    xyarc = Graphics3D[{blue, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 0, 
                \(polar[XY]\)[[2
                  ]],  .1}]]}]; 
Clear[xyarc1]; 

    xyarc1 = Graphics3D[{color1, 
          Line[Table[Cos[t]\ X0 + Sin[t]\ Y0, {t, \(- .5\), 0,  .1}]]}]; 

    Clear[xyarc2]; 

    xyarc2 = Graphics3D[{color1, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 
                \(polar[XY]\)[[2]], 
                \(polar[XY]\)[[2]] + 
                  0.3,  .1}]]}]; 
Clear[stick]; 

    stick = Graphics3D[{red, Line[{scale\ Z, XY}]}]; 

    \(xaxislabel = Text["x", {0.363, 0.26}]; \)\)yaxislabel = Text["y", {0.73, 0.3}]; zaxislabel = Text["z", { .5, 0.67}]; stickplot = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; rho = Text["r", {0.53, 0.45}]; rlabel = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; phi = Text["s", {0.55, 0.57}]; slabel = 
      Show[xyplane, zpoint, axes, ray, point, zarc, zarc1, zarc2, 
        Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {phi, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; theta = Text["t", { .5, 0.3}]; tlabel = 
      Show[xyplane, axes, xypoint, xpoint, ray, stick, point, xyarc, xyarc1, 
        xyarc2, Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; all = 
      Show[xyplane, axes, ray, stick, point, zpoint, xpoint, xypoint, xyarc, 
        xyarc1, xyarc2, zarc, zarc1, zarc2, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, phi, rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity];

Crisscross

\(Clear[x]; 

    crisscross = 
      Plot[{Cos[x], Sin[x]}, {x, π\/8, \(3\ π\)\/8}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.02]}}, 
        Axes -> None, AspectRatio -> 1\/3, 
        DisplayFunction -> Identity]; \)

B.1)  Sources, sinks, and Gauss's formula in 3D

B.1.a)

Given a three-dimensional vector field 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]}:

\(Clear[Field, m, n, p, x, y, z]; 

    Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

You calculate the divergence, divField[x, y, z], of Field[x, y, z] via the formula:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]\)

Not much of surprise. 
What do you do with this formula?

You use it the same way you use its 2D counterpart.

B.1.b.i) Gauss's 3D formula

Gauss's 3D formula says that if you go with a 3D vector field 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]}, 
and if R is a solid body in three dimensions with boundary surface (skin) C, then
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\ 
            divField[x, y] dx dy 
              dz\)\)\)
        = 
        (∫_^)(∫_( C)) 
            Field . outerunitnormal dA
        = flow of Field[x, y, z] across C.
Use Gauss's formula to measure the flow of the vector field 
     Field[x, y, z] = {y -  0.5, x y, (z^2)}
across the surface of the box consisting of all points {x, y, z} with
       -1 ≤ x ≤ 3, 0 ≤ y ≤ 2, and -2 ≤ z ≤ 4.

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {y - 0.5, x\ y, z\^2}; 
    
Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Gauss's formula tells you that the net flow of Field[x, y, z] across the surface of the box with
       -1 ≤ x ≤ 3, 0 ≤ y ≤ 2, and -2 ≤ z ≤ 4
is 
       (∫_(-2)^4)
          (∫_0^2)
            (∫_(-1)^3) divField[x, y, z] dz 
              dy dx
Here you go:

\(Clear[divF]; 

    divField[x_, y_, z_] = 
      D[m[x, y, z], x] + \ D[n[x, y, z], y] + D[p[x, y, z], z]; 

    ∫\_\(-2\)\%4
        \(∫\_0\%2
          \(∫\_\(-1\)\%3 divField[x, y, z] \[DifferentialD]z 
            \[DifferentialD]y \[DifferentialD]x\)\)\)

Very positive. 
Strong flow from inside the box to outside the box. 
There must be lots of spigots turned on inside this box.

B.1.b.ii) Sources and sinks

In three dimensions, if 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]},
then
       \(divField[x, y, z] = 
	\ \ 
        D[m[x, y, z], x] + \ D[n[x, y, z], y] + \ D[p[x, y, z], z]\) 
tells you about sources and sinks the same way that divField[x, y] tells you about sources and sinks in two dimensions: 
-> The points {x, y, z} with divField[x, y, z] > 0 are sources of new fluid.  
-> The points {x, y, z} with divField[x, y, z] < 0 are sinks for old fluid.
Use Gauss's formula to explain why this interpretation is legitimate.

Here's how to see why:
Take a small sphere C centered at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}.  Calculate the flow across C which, according to Gauss's formula, is 
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\ 
            divField[x, y] dx dy 
              dz\)\)\)
where R is the solid region consisting of the sphere C and everything inside C.
Here's the kicker:  
If you start out with
       divField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] > 0, 
then divField[x, y, z] is positive for all {x, y, z}'s close to {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}, so if C is small enough, 
       divField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] > 0
for all {x, y, z}'s inside C.  This tells you that that for small spheres C centered at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])},
       \(flow - across - C = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\)\ 
                divField[x, y] dx dy 
                  dz\)\) > 0\).
And this tells you that if divField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] > 0, then the net flow of Field[x, y, z] across small spheres centered at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])} is from inside to outside.  
The upshot:
If divField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] > 0, then the point{(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])} is a source of new fluid.
Similarly, if divField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] < 0, then the net flow of Field[x, y, z] across small spheres centered at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])} is from outside to inside.  Consequently, if divField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] < 0, then the point {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])} is a sink (or drain) for old fluid.

B.1.b.iii)

Given 
       Field[x, y, z] = {y -  0.5, x y, (z^2)}, 
say how to identify the points {x, y, z} that are sources, and the points {x, y, z} that are sinks.

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {y - 0.5, x\ y, z\^2}
        \)Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]};

Calculate the divergence:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]\)

Find out where divField[x, y, z] is 0:

Solve[divField[x, y, z] == 0]

This says that the points {x, y, z} with x = -2 z are neither sources nor sinks.
Now look at:

\(Clear[a]; 
divField[x, y, z] /. x -> \(-2\)\ z + a\)

Food for thought. 
Upon reflection, this tells you that if 
       x > -2 z, then divField[x, y, z] > 0. 
and if 
       x < -2 z, then divField[x, y, z] < 0.
As a result, the points {x, y, z} that are sources are the points with x > -2 z, and the points {x, y, z} that are sinks are the points with x < -2 z.

B.1.b.iv)

Given 
       Field[x, y, z] = {x -  z, y -  x, z -  y},  
say how to identify the points {x, y, z} that are sources, and the points {x, y, z} that are sinks.
Use your answer to determine whether the net flow of this vector field across the sphere of radius 2 centered at {4, 2, 1} is from inside to outside or from outside to inside.

\(Clear[m, n, p, Field, x, y, z]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {x - z, y - x, z - y}; 
    
Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

divField[x, y, z] is given by:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      ∂\_x m[x, y, z] + ∂\_y n[x, y, z] + 
        ∂\_z p[x, y, z]\)

No matter what {x, y, z} happens to be,
       divField[x, y, z] = 3.
So:
All points {x, y, z} are sources.
And now without further calculation, you know that the net flow of this vector field across sphere of radius 2 centered at {4, 2, 1} is from inside to outside.
Reason: There are no sinks inside the sphere to absorb excess outside-to-inside flow.

B.1.c)

Explain the reasoning behind Gauss's formula

A heavy, long-winded explanation is possible, but leave it at this:
Once you have a good feeling for the ideas in B.2), B.3), and B.4), you can build your own explanation of Gauss's formula by copying and pasting the explanation of the two dimensional Gauss-Green formula from an earlier lesson and making technical adjustments. 
If you cannot find satisfaction without seeing some more details, see:

B.2)  Measuring area on surfaces

B.2.a) Area of a parallelogram in three dimensions

Here's a parallelogram in 3 dimensions:

basepoint = {1, 2, 0}; X = {-2, 1, 1}; Y = {1, -2, 1}; parallelogram = 
      Show[Graphics3D[
          Polygon[{basepoint, basepoint + X, basepoint + X + Y, 
              basepoint + Y}]], ViewPoint -> CMView, Axes -> True, 
        AxesLabel -> {"x", "y", "z"}];

Because of its unfortunate position, it might seem hard to measure the area of this parallelogram, but the cross product can bail you out.  
How?

Throw in the vectors that define this parallelogram:

Show[parallelogram, 
      Arrow[X, Tail -> basepoint, VectorColor -> Blue], 
      Arrow[Y, Tail -> basepoint, VectorColor -> Blue], 
      Graphics3D[{PointSize[0.03], Point[basepoint]}], 
      Graphics3D[Text["X", basepoint + (X/2) + {0, 0.2, 0}]], 
      Graphics3D[Text["Y", basepoint + (Y/2) + {0.2, 0, 0}]], 
      ViewPoint -> CMView, PlotRange -> All, Axes -> True, 
      AxesLabel -> {"x", "y", "z"}];

As you probably remember,
       ‖X×Y‖ = 
        ‖X‖ 
          ‖Y‖ 
          |Sin[angle between]|.
This quantity is also equal to the area of this parallelogram; so the area in square units of this unfortunately positioned parallelogram is given by:

XcrossY = X\[Cross]Y; (Sqrt[XcrossY . XcrossY])

Easy.
Rerun for other choices of X and Y to get the hang of it.

B.2.b.i)  Measuring area on planes

Here's a look at part of the plane
     x +  y +  z = 2:

\(Clear[x, y, z, u, v]; 
x[u_, v_] = u; 
\(y[u_, v_] = v; \)\)z[u_, v_] = 2 - u - v; surface = 
      ParametricPlot3D[{x[u, v], y[u, v], z[u, v]}, {u, 0, 1}, {v, 0, 1}, 
        PlotPoints -> {2, 2}, Axes -> Automatic, 
        AxesLabel -> {"x", "y", "z"}, 
        ViewPoint -> CMView];

Here is another way of getting a plot of the same surface:
This plot also includes the vectors that generate the parallelogram.

{a, b} = {0, 0}; h = 1; Show[
      Graphics3D[
        Polygon[{{x[a, b], y[a, b], z[a, b]}, {x[a + h, b], y[a + h, b], 
              z[a + h, b]}, {x[a + h, b + h], y[a + h, b + h], 
              z[a + h, b + h]}, {x[a, b + h], y[a, b + h], z[a, b + h]}}]], 
      Arrow[{x[a + h, b], y[a + h, b], z[a + h, b]} - {x[a, b], y[a, b], 
            z[a, b]}, Tail -> {x[a, b], y[a, b], z[a, b]}, 
        VectorColor -> Blue], 
      Arrow[{x[a, b + h], y[a, b + h], z[a, b + h]} - {x[a, b], y[a, b], 
            z[a, b]}, Tail -> {x[a, b], y[a, b], z[a, b]}, 
        VectorColor -> Blue], PlotRange -> All, 
      Axes -> Automatic, AxesLabel -> {"x", "y", "z"}, 
      ViewPoint -> CMView];

For this choice of a, b, and h, the corresponding parallelogram has area measurement:

X = {x[a + h, b], y[a + h, b], z[a + h, b]} - {x[a, b], y[a, b], 
          z[a, b]}; Y = {x[a, b + h], y[a, b + h], z[a, b + h]} - {x[a, b], y[a, b], 
          z[a, b]}; XcrossY = X\[Cross]Y; (Sqrt[XcrossY . XcrossY])

For unspecified choices of a, b, and h, the corresponding parallelogram has area measured by:

\(Clear[a, b, h]; 

    X = {x[a + h, b], y[a + h, b], z[a + h, b]} - {x[a, b], y[a, b], 
          z[a, b]}; 

    \(Y = {x[a, b + h], y[a, b + h], z[a, b + h]} - {x[a, b], y[a, b], 
          z[a, b]}; \)\)XcrossY = X\[Cross]Y; (Sqrt[XcrossY . XcrossY])

What is the area conversion factor 
       (Subscript[SA, xyz])[u, v]
that you use to convert uv-paper area measurement into xyz-surface area measured on the plane x +  y +  z = 2?

The uv-paper rectangle with corners at 
       {a, b}, {a +  h, b}, {a +  h, b +  h}, and {a, b +  h} 
has uv-paper area given by:

\(Clear[h]; 
uvarea = h\^2\)

This uv-paper rectangle plots out in xyz-coordinates as the parallelogram with corners at 
       {x[a, b], y[a, b], z[a, b]},
       {x[a +  h, b], y[a +  h, b], z[a +  h, b]},
       {x[a +  h, b +  h], y[a +  h, b +  h], z[a +  h, b +  h]}, 
and
       {x[a, b +  h], y[a, b +  h], z[a, b +  h]}.
The area of this parallelogram is given by:

\(Clear[a, b, h]; 

    X = {x[a + h, b], y[a + h, b], z[a + h, b]} - {x[a, b], y[a, b], 
          z[a, b]}; 

    \(Y = {x[a, b + h], y[a, b + h], z[a, b + h]} - {x[a, b], y[a, b], 
          z[a, b]}; \)\)XcrossY = X\[Cross]Y; planararea = (Sqrt[XcrossY . XcrossY])

The area conversion factor (Subscript[SA, xyz])[u, v] that you integrate to convert uv-paper area into xyz-surface area measured on the plane x +  y +  z = 2 is:

\(Clear[SAxyz]; 
SAxyz[u_, v_] = PowerExpand[planararea\/uvarea]\)

Still pretty easy.

B.2.b.ii)  Measuring planar area

Here is a plot of the uv-paper circle 
       (u^2) +  (v^2) = 4  
on the plane 
       x +  y +  z = 2
under the parameterization
       x[u, v] = u,
       y[u, v] = v, and
       z[u, v] = 2 -  u -  v 
used above.

\(Clear[x, y, z, u, v, t]; 
x[u_, v_] = u; 
\(y[u_, v_] = v; \)\)z[u_, v_] = 2 - u - v; u[t_] = 2 Cos[t]; v[t_] = 2 Sin[t]; curve = 
      ParametricPlot3D[{x[u[t], v[t]], y[u[t], v[t]], z[u[t], v[t]]}, {t, 0, 
          2 π}, DisplayFunction -> Identity]; {a, b} = {-4, -4}; h = 8; \(plane = 
      Graphics3D[
        Polygon[{{x[a, b], y[a, b], z[a, b]}, {x[a + h, b], y[a + h, b], 
              z[a + h, b]}, {x[a + h, b + h], y[a + h, b + h], 
              z[a + h, b + h]}, {x[a, b + h], y[a, b + h], z[a, b + h]}}]]; 

    
Show[curve, plane, Axes -> Automatic, 
      AxesLabel -> {"x", "y", "z"}, 
      ViewPoint -> CMView, DisplayFunction -> $DisplayFunction]; 
    \)

Come up with a measurement of area enclosed by this curve as measured on the plane
       x +  y +  z = 2.

Area on the plane is (Subscript[SA, xyz])[u, v] times uv-paper area as calculated above in part i):

SAxyz[u, v]

The circle (u^2) +  (v^2) = 4 is a circle of radius 2 on uv-paper. It encloses an area of (2^2) π = 4  π square units measured on uv-paper. 
Its plot on the plane x +  y +  z = 2 encloses a total of:

SAxyz[u, v] 4 π

in square units.

B.2.c.i)  Measuring area on curved surfaces

Here's a portion of the surface whose parametric equations are:

\(Clear[x, y, z, u, v]; 
x[u_, v_] = u; 
\(y[u_, v_] = v; \)\)z[u_, v_] = u v; {ulow, uhigh} = {0, 1}; {vlow, vhigh} = {0, 2}; surface = 
      ParametricPlot3D[
        Evaluate[{x[u, v], y[u, v], z[u, v]}], {u, ulow, uhigh}, {v, vlow, 
          vhigh}, ViewPoint -> CMView, BoxRatios -> Automatic, 
        Axes -> Automatic, 
        AxesLabel -> {"x", "y", "z"}];

Measure the surface area of the plotted portion of this surface.

Mathematica plots this surface by approximating the true surface with a whole herd of little patches that resemble parallelograms.
The message is clear.
You can use what you know about measuring the area of parallelograms to help to measure the area of curved surfaces.  
See what some parallelograms will do:

\(Clear[h, basepoint, leg1, leg2, parallelogram]; 

    basepoint[u_, v_] = {x[u, v], y[u, v], z[u, v]}; 

    \(leg1[u_, v_, h_] = {x[u + h, v], y[u + h, v], z[u + h, v]} - 
        basepoint[u, v]; \)\)leg2[u_, v_, h_] = {x[u, v + h], y[u, v + h], z[u, v + h]} - 
        basepoint[u, v]; parallelogram[u_, v_, h_] = 
      Graphics3D[{Black, 
          Polygon[{basepoint[u, v], basepoint[u, v] + leg1[u, v, h], 
              basepoint[u, v] + leg1[u, v, h] + leg2[u, v, h], 
              basepoint[u, v] + leg2[u, v, h]}]}]; jump = 0.3; h = 0.25; \(parallelograms = 
      Table[parallelogram[u, v, h], {u, ulow, uhigh - jump, jump}, {v, vlow, 
          vhigh - jump, jump}]; 


    Show[surface, parallelograms, ViewPoint -> CMView, 
      Axes -> Automatic, AxesLabel -> {"x", "y", "z"}]; 
    \)

Hug that surface, mama! 
The smart money says:
The smaller the parallelograms are, the closer they hug the surface.
Check it out:

jump = 0.2; h = 0.1; \(parallelograms = 
      Table[parallelogram[u, v, h], {u, ulow, uhigh - jump, jump}, {v, vlow, 
          vhigh - jump, jump}]; 


    Show[surface, parallelograms, ViewPoint -> CMView, 
      Axes -> Automatic, AxesLabel -> {"x", "y", "z"}]; \)

The smart money wins. 
These smaller parallelograms are really glued to the surface.  Even smaller parallelograms will share even more ink with the surface. 
Here's how to get the area conversion factor (Subscript[SA, xyz])[u, v] for going from uv-paper area to xyz-area on the surface at a uv-point {u, v}: 
Notice that for small positive h, the area conversion factors at a uv-point {u, v} are nearly the same for area on the surface and for area on the parallelograms plotted above which are determined by the vectors:

\(Clear[h]; 
leg1[u, v, h]\)leg2[u, v, h]

For fixed {u, v} and h, the area conversion factor for converting uv-area to xyz-area on the parallelogram determined by leg1[u, v, h] and leg2[u, v, h] is given by:

cross = leg1[u, v, h]\[Cross]leg2[u, v, h]; \(Clear[paraareafactor]; 

    paraareafactor[u_, v_, h_] = \@\(cross . cross\)\/h\^2\)

The area conversion factor (Subscript[SA, xyz])[u, v] for the surface at {u, v} is the limiting case of the above as h -> 0:

\(Clear[SAxyz]; 

    SAxyz[u_, v_] = Limit[paraareafactor[u, v, h], h -> 0]\)

The total surface area (in square units) of the surface plotted above is
       (∫_vlow^vhigh)(∫_ulow^uhigh) 
          (Subscript[SA, xyz])[u, v] du dv:

surfarea = 
      NIntegrate[SAxyz[u, v], {u, ulow, uhigh}, {v, vlow, vhigh}]

And once you've done one of these, you've done them all.

B.2.c.ii)

Here is a direct formula for the area conversion factor 
(Subscript[SA, xyz])[u, v] that converts uv-paper area measurements into area measurements on the surface whose parametric formulas are {x[u, v], y[u, v], z[u, v]}:

\(Clear[x, y, z, u, v, cross, SAxyz]; 

    cross[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]; 
SAxyz[u_, v_] = \@\(cross[u, v] . cross[u, v]\)\)

Nasty looking, but useful.
Explain where this formula for (Subscript[SA, xyz])[u, v] comes from.

You can make a good run at the formula by clearing the parametric functions x[u, v], y[u, v], and z[u, v], and doing what was done in the last part above:

\(Clear[x, y, z, u, v, h, basepoint, leg1, leg2, parallelogram]; 

    basepoint[u_, v_] = {x[u, v], y[u, v], z[u, v]}\)

leg1[u_, v_, h_] = {x[u + h, v], y[u + h, v], z[u + h, v]} - 
        basepoint[u, v]

leg2[u_, v_, h_] = {x[u, v + h], y[u, v + h], z[u, v + h]} - 
        basepoint[u, v]

As h closes in on 0, the area conversion factor for the parallelogram determined by the vectors:

leg1[u, v, h]

and

leg2[u, v, h]

with tails at

basepoint[u, v]

closes in on the area conversion factor (Subscript[SA, xyz])[u, v] for the surface at {u, v}.
The area conversion factor for this parallelogram is given by
       ‖leg1[u, v, h]×leg2[u, v, h]
                  (‖/(h^2)).
As a result, the formula for the area conversion factor surface at {u, v} is
       (lim_(h  ->  0))‖leg1[u, v, h]×leg2[u, v, h]
                    (‖/(h^2)).
This is the same as
       (lim_(h  ->  0))‖leg1[u, v, h]/h
                  ×leg2[u, v, h]/h
                  ‖.
This is the same as
       ‖(lim_(h  ->  0))(leg1[u, v, h]/h)×(lim_(h  ->  0))(((leg2[u, v, h])/h))‖.
Look at:

leg1[u, v, h]/h

leg2[u, v, h]/h

(lim_(h  ->  0)) (leg1[u, v, h]/h) is given by:

first = {D[x[u, v], u], D[y[u, v], u], D[z[u, v], u]}

(lim_(h  ->  0)) (leg2[u, v, h]/h) is given by:

second = {D[x[u, v], v], D[y[u, v], v], D[z[u, v], v]}

The area conversion factor (Subscript[SA, xyz])[u, v] for the surface is the length of cross product:

\(Clear[cross, SAxyz]; 
cross[u_, v_] = first\[Cross]second; 

    SAxyz[u_, v_] = \@\(cross[u, v] . cross[u, v]\)\)

The final formula is a mess to remember, but you can get the same formula quickly via:

\(Clear[cross, SAxyz]; 

    cross[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]; 
SAxyz[u_, v_] = \@\(cross[u, v] . cross[u, v]\)\)

And this is the explanation.

B.2.c.i)

Measure the surface area of the surface of a sphere of radius r.

Go to spherical coordinates.

Show[all, DisplayFunction -> $DisplayFunction];

The surface of the sphere of radius r is described by 
       \({x[s, t], y[s, t], z[s, t]} = 
	\ \ \ \ \ 
        \ {r\ Sin[s] Cos[t], r\ Sin[s] Sin[t], r\ Cos[s]}\) 
with s running from 0 to π and t running from 0 to 2  π.

\(Clear[x, y, z, r, s, t]; 
    
{x[s_, t_], y[s_, t_], z[s_, t_]} = {r\ Sin[s]\ Cos[t], 
        r\ Sin[s]\ Sin[t], r\ Cos[s]}\)

To measure the area of surface of the sphere of radius r, all you have to do is calculate 
       (∫_0^π)(∫_0^(2  π)) 
          (Subscript[SA, xyz])[s, t] dt ds:

\(Clear[SAxyz, cross]; 

    cross[s_, t_] = 
      D[{x[s, t], y[s, t], z[s, t]}, s]\[Cross]D[{x[s, t], y[s, t], z[s, t]}, 
          t]; 
SAxyz[s_, t_] = 
      √TrigExpand[cross[s, t] . cross[s, t]]\)

(∫_0
        ^π)(∫_0^(2 π))SAxyz[s, t] \[DifferentialD]t 
          \[DifferentialD]s

If your version of Mathematica gave you 0, then your version of Mathematica screwed up. Here is another way of getting the correct answer:

2 (∫_0
          ^((π/2)))(∫_0^(2 π))SAxyz[s, t] 
            \[DifferentialD]t \[DifferentialD]s

4  π (r^2) is correct.

B.3)  Surface integrals

Just as you can integrate functions on curves with respect to arc length, you can integrate functions on surfaces with respect to surface area.
And you don't have to go to a lot of trouble to do it.

B.3.a)

Here is a parameterized surface C:

\(Clear[x, y, z, u, v]; 
x[u_, v_] = u\ Cos[v]; 

    \(y[u_, v_] = \(-u\)\ Sin[v]; \)\)z[u_, v_] = v; {{ulow, uhigh}, {vlow, vhigh}} = {{1, 3}, {0, ((3 π)/2)}}; 
    surface = 
      ParametricPlot3D[{x[u, v], y[u, v], z[u, v]}, {u, ulow, uhigh}, {v, 
          vlow, vhigh}, ViewPoint -> CMView, Boxed -> False, 
        AxesLabel -> {"x", "y", "z"}];

Calculate 
       (∫_^)(∫_( R)^) 
          ((x^2) +  (y^2)) dA
where the integral is taken with respect to surface area.

(∫_^)(∫_( R)^) 
          ((x^2) +  (y^2)) dA
is calculated via the formula
       (∫_vlow^vhigh)(∫_ulow^uhigh) 
          (x[u, v]^2 +  y[u, v]^2) (Subscript[SA, xyz])[u, v] du 
            dv.
Calculating this is duck soup.
First you need the area conversion factor (Subscript[SA, xyz])[u, v]:

\(Clear[SAxyz]; 

    cross = D[{x[u, v], y[u, v], z[u, v]}, u]
          \[Cross]D[{x[u, v], y[u, v], z[u, v]}, v]; 

    SAxyz[u_, v_] = \@TrigExpand[cross . cross]\)

Now turn the integral
       (∫_vlow^vhigh)(∫_ulow^uhigh) 
          (x[u, v]^2 +  y[u, v]^2) (Subscript[SA, xyz])[u, v] du 
            dv.
over to the machine:

integral = 
      NIntegrate[
        (x[u, v]^2 + y[u, v]^2) SAxyz[u, v], {u, ulow, uhigh}, {v, 
          vlow, vhigh}]

Not exciting, but not hard.

B.4)  Surface integrals for measuring flow across surfaces

B.4.a)

Here's a surface:

\(Clear[x, y, z, u, v]; 
x[u_, v_] = u; 
\(y[u_, v_] = v; \)\)z[u_, v_] = (1/3) (10 - (u^2) - (v^2)); {ulow, uhigh} = {-2, 2}; {vlow, vhigh} = {-2, 2}; surface = 
      ParametricPlot3D[
        Evaluate[{x[u, v], y[u, v], z[u, v]}], {u, ulow, uhigh}, {v, vlow, 
          vhigh}, AxesLabel -> {"x", "y", "z"}, 
        ViewPoint -> CMView];

Take {a, b} = {0.5, 1.0} and look at this display of the two curves
       (Subscript[P, 1])[u] = {x[u, b], y[u, b], z[u, b]}
with ulow ≤ u ≤ uhigh
and
       (Subscript[P, 2])[v] = {x[a, v], y[a, v], z[a, v]} 
with vlow ≤ v ≤ vhigh
together with the surface plotted above:

{a, b} = {0.5, 1.0}; \(Clear[P1, P2]; 
P1[u_] = {x[u, b], y[u, b], z[u, b]}; 

    \(P2[v_] = {x[a, v], y[a, v], z[a, v]}; \)\)curve1 = 
      ParametricPlot3D[Evaluate[P1[u]], {u, ulow, uhigh}, 
        DisplayFunction -> Identity]; curve2 = 
      ParametricPlot3D[Evaluate[P2[v]], {v, vlow, vhigh}, 
        DisplayFunction -> Identity]; \(surface = Insert[surface, EdgeForm[], {1, 1}]; 


    curves = Show[surface, curve1, curve2, Boxed -> False, 
        ViewPoint -> CMView, BoxRatios -> Automatic, 
        DisplayFunction -> $DisplayFunction]; \)

Two intersecting curves running on the surface.
Now calculate the area conversion factor (Subscript[SA, xyz])[u, v]

\(Clear[SAxyz, cross]; 

    cross[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]; 
SAxyz[u_, v_] = \@\(cross[u, v] . cross[u, v]\)\)

Keeping the same {a, b} as above, add to the last plot the three vectors
       D[{x[u, v], y[u, v], z[u, v]}, u] /. {u -> a, v -> b},
       D[{x[u, v], y[u, v], z[u, v]}, v]}] /. {u -> a, 
          v -> b},and 
       cross[a, b],
all with their tails at the point of intersection:

tail = {x[a, b], y[a, b], z[a, b]}; \(vector1 = 
      D[{x[u, v], y[u, v], z[u, v]}, u] /. {u -> a, v -> b}; 

    vector2 = 
      D[{x[u, v], y[u, v], z[u, v]}, v] /. {u -> a, v -> b}; 

    vector3 = cross[a, b]; 


    Show[curves, Arrow[vector1, Tail -> tail], 
      Arrow[vector2, Tail -> tail], Arrow[vector3, Tail -> tail]]; 
    \)

Totally radical or what!
Explain why this happened and how it relates to the area conversion factor (Subscript[SA, xyz])[u, v].

Both
       D[{x[u, v], y[u, v], z[u, v]}, u], and     
       D[{x[u, v], y[u, v], z[u, v]}, v]
are tangent to the surface at {x[u, v], y[u, v], z[u, v]},  
so their cross product, cross[u, v], is automatically normal (perpendicular) to the surface at {x[u, v], y[u, v], z[u, v]}. 
The upshot:
The area conversion factor 
       (Subscript[SA, xyz])[u, v] = (Sqrt[cross[u, v] . cross[u, v]])
is the same as the length of the normal vector 
       cross[u, v].

B.4.b)

Given a 3D vector field 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]},
and given a surface parameterized by 
     {x[u, v], y[u, v], z[u, v]} 
with ulow ≤ u ≤ uhigh and vlow ≤ v ≤ vhigh, 
you can measure the flow of Field[x, y, z] across the surface by calculating
       (∫_vlow^vhigh)(∫_ulow^uhigh) 
          Field[x[u, v], y[u, v],  z[u, v]] . normal[u, v] du 
            dv
where normal[u, v] = cross[u, v] is given by:

\(Clear[x, y, z, u, v, normal]; 

    normal[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]\)

Where does the flow-across formula  
       (∫_vlow^vhigh)
          (∫_ulow
            ^uhigh) Field[x[u, v], y[u, v],  z[u, v]] . normal[u, v] 
            du dv
come from?

Put normcompField[x[u, v], y[u, v], z[u, v]] equal to the component of Field[x[u, v], y[u, v], z[u, v]] in the direction perpendicular to the surface at {x[u, v], y[u, v], z[u, v]}.
You measure the flow of Field[x, y, z] across the surface by calculating
       flowacross = 
        (∫_^)(∫_( surface)^) 
            normcompField[x[u, v], y[u, v], z[u, v]] dA  
        = 
        (∫_vlow^vhigh)(∫_ulow^uhigh) 
            normcompField[x[u, v], y[u, v],  z[u, v]] (Subscript[SA, xyz])[u, v] 
              du dv.
Put this in your pocket for a minute. 
Remember that normal[u, v] as calculated above is perpendicular to the surface at the point {x[u, v], y[u, v], z[u, v]}. At this point, the component of Field[x[u, v], y[u, v], z[u, v]] in the direction perpendicular to the surface is
       normcompField[x[u, v], y[u, v], z[u, v]]
       =Field[x[u, v], y[u, v], z[u, v]].normal[u, v]/(Sqrt[normal[u, v] . normal[u, v]]).
       
Because (Subscript[SA, xyz])[u, v] = (Sqrt[normal[u, v] . normal[u, v]]),
       flowacross = 
        (∫_^)(∫_( surface)^) 
            normcompField[x[u, v], y[u, v], z[u, v]] dA  
       
        = 
        (∫_vlow^vhigh)(∫_ulow^uhigh) 
            normcompField[x[u, v], y[u, v],  z[u, v]] (Subscript[SA, xyz])[u, v] 
              du dv
              =(∫_vlow^vhigh)(∫_ulow^uhigh) Field[x[u, v], y[u, v],  z[u, v]].(normal[u, v]/(Subscript[SA, xyz])[u, v])(Subscript[SA, xyz])[u, v]dudv
              
       =(∫_vlow^vhigh)(∫_ulow^uhigh) Field[x[u, v], y[u, v],  z[u, v]].(normal[u, v]/(Subscript[SA, xyz])[u, v])dudv
       
        = 
        (∫_vlow^vhigh)
            (∫_ulow
              ^uhigh) Field[x[u, v], y[u, v],  z[u, v]] . normal[u, v] 
              du dv.
because the (Subscript[SA, xyz])[u, v] terms miraculously cancel out.
Explanation finished.

B.4.c)

Here's a new surface:

\(Clear[x, y, z, u, v]; 
x[u_, v_] = v; 

    \(y[u_, v_] = 2\ E\^\(-v\)\ Sin[π\ u]; \)\)z[u_, v_] = u; {ulow, uhigh} = {0, 2}; {vlow, vhigh} = {0.5, 3.5}; surface = 
      ParametricPlot3D[
        Evaluate[{x[u, v], y[u, v], z[u, v]}], {u, ulow, uhigh}, {v, vlow, 
          vhigh}, ViewPoint -> CMView, 
        AxesLabel -> {"x", " y", "z"}];

Here is this surface together with some of the normal vectors whose lengths measure the area conversion factor 
       (Subscript[SA, xyz])[u, v], 
which converts uv-paper area measurements into xyz-area measurements on the surface:

\(Clear[normal]; 

    normal[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]; 
normals = 
      Table[Arrow[normal[u, v], Tail -> {x[u, v], y[u, v], z[u, v]}], {
          u, 0.5, 1.5, 0.5}, {v, 1.5, 2.5}]; 

    surface = Insert[surface, EdgeForm[], {1, 1}]; 


    Show[surface, normals, ViewPoint -> CMView, Boxed -> False, 
      AxesLabel -> {"x", "y", "z"}]; \)

The normals are pointing out from the side of the surface you are sitting on.
Determine whether the net flow of the vector field
       Field[x, y, z] = {y z, x z, x y} 
across this surface is with the plotted normals or is against the plotted normals.

Enter the vector field:

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {y\ z, x\ z, x\ y}; 

    Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Reenter the parameterization of the surface, and calculate normal[u, v]:

\(Clear[u, v]; 
x[u_, v_] = v; 

    \(y[u_, v_] = 2\ E\^\(-v\)\ Sin[π\ u]; \)\)z[u_, v_] = u; {ulow, uhigh} = {0, 2}; {vlow, vhigh} = {0.5, 3.5}; \(Clear[normal]; 

    normal[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]\)

Calculate  
       (∫_vlow^vhigh)
          (∫_ulow^uhigh) Field[x[u, v], y[u, v]] . normal[u, v] 
            du dv:

NIntegrate[
      Field[x[u, v], y[u, v], z[u, v]] . normal[u, v], {u, ulow, uhigh}, {v, 
        vlow, vhigh}, AccuracyGoal -> 2]

Positive.
This tells you that the net flow of Field[x, y, z] across this surface is in the direction of the plotted normals.
Take another look:

Show[surface, normals, ViewPoint -> CMView, Boxed -> False, 
      AxesLabel -> {"x", "y", "z"}];

In short, the net flow of 
       Field[x, y, z] = {y z, x z, x y} 
across this surface is from behind the surface to the front of the surface.

B.5)  Flux of the electric field

Way, way back, at the end of the (18^th) century, a fellow named Coulomb played with electrostatic problems.  He rubbed cats with amber rods, and had an all-around good time with his science.  More than that, he managed to do some calculations of forces that electrically charged objects exert on each other.
  
He found out that in free open space, the force that a particle located at {x, y, z} with charge q[x, y, z] exerts on a particle located at {xx, yy, zz} with charge q[xx, yy, zz] is given by:

\(Clear[K, q, x, y, z, xx, yy, zz]; 
R = {x, y, z} - {xx, yy, zz}; 

    Force = \(K\ q[x, y, z]\ q[xx, yy, zz]\ R\)\/\((R . R)\)\^\(3/2\)\)

Note that R runs parallel to the vector that points from {xx, yy, zz} to {x, y, z}.

Look at the magnitude of this force:

(Sqrt[Together])[Force . Force]

This is the same as:

((K q[x, y, z] q[xx, yy, zz])/Expand)[R . R]

This means that the force is directed on the line between the two points, and its magnitude is 
        ((K q[x, y, z] 
              q[xx, yy, zz])/(((distance between the two points))^2))
              
because the denominator
       R . R = (((distance between the two points))^2).

That constant K involves several other physical quantities and ultimately depends only on the system of units used in the calculation.  As usual in mathematics, you specify standard math course units of measure by adjusting the units so that K = 1, and then you leave it to your physics friends to tell you about the K they actually use.    
For visualizations of Coulomb's law, electrical folks like the idea of the electric field. The three dimensional electric field of assumes a charge 
       q[x, y, z] = 1 at {x, y, z}, 
and thus is the force per unit charge at {x, y, z}:

\(Clear[ElectField]; 

    ElectField[x_, y_, z_] = Force /. {K -> 1, q[x, y, z] -> 1}\)

Here is a plot of the 3D electric field for {xx, yy, zz} = {0, 0, 0} and q[xx, yy, zz] = 6:

{xx, yy, zz} = {0, 0, 0}; q[xx, yy, zz] = 6; Show[Graphics3D[{Red, PointSize[0.06], Point[{xx, yy, zz}]}], 
      Table[Arrow[ElectField[x, y, z], Tail -> {x, y, z}, 
          VectorColor -> Red], {x, -3, 3, 2}, {y, -3, 3, 2}, {z, 
          -3, 3, 2}]];

The field vectors point in the direction of greatest voltage drop because the electric field is the gradient of the negative voltage.

B.5.a)

Measure the flux across the surface of the sphere 
       (x^2) + (y^2) + (z^2) = 1
resulting from a single charge of strength 3 placed at {0, 0, 0}.
Then take any positive radius r and calculate the flux across the surface of the sphere
       (x^2) + (y^2) + (z^2) = (r^2)
of the electric field plotted above.
Discuss any noteworthy outcomes.

The flux of the electric field across a surface is just another name for the flow of the electric field across that surface.
Load in the electric field:

Clear[K, q, x, y, z, a]; R = {x, y, z} - {0, 0, 0}; Force = ((K q[x, y, z] q[0, 0, 0] R)/(((R . R))^(3/2)))

q[0, 0, 0] = 3; \(Clear[ElectField]; 

    ElectField[x_, y_, z_] = 
      Force /. {K -> 1, q[x, y, z] -> 1}\)

-> To measure the flux through (x^2) + (y^2) + (z^2) = 1:
Use spherical coordinates for the sphere (x^2) + (y^2) + (z^2) = 1:

\(Clear[x, y, z, r, s, t, normal]; 
r = 1; 

    \({x[s_, t_], y[s_, t_], z[s_, t_]} = {r\ Sin[s]\ Cos[t], 
        r\ Sin[s]\ Sin[t], r\ Cos[s]}; \)\){{slow, shigh}, {tlow, thigh}} = {{0, π}, {0, 2 π}}; normal[s_, t_] = 
      TrigExpand[
        D[{x[s, t], y[s, t], z[s, t]}, s]
            \[Cross]D[{x[s, t], y[s, t], z[s, t]}, t]]

You can see that these normals point out from the sphere because the third slot is positive for 0 < s < (π/2) and the third slot is negative for (π/2) < s < π. 
Confirm with a plot:

Show[
      Table[Arrow[normal[s, t], Tail -> {x[s, t], y[s, t], z[s, t]}, 
          VectorColor -> Blue], {s, (π/4), ((3 π)/4), 
          (π/4)}, {t, (π/4), ((7 π)/4), (π/2)}], 
      ViewPoint -> CMView, Axes -> True, 
      AxesLabel -> {"x", "y", "z"}];

Good.
The flux of an electric field across a surface is just the flow of the same electric field across the surface.
See B.4) to learn why this is given by 
       (∫_^)(∫_( surface)^) 
          ElectField . unitnormal dA
        = 
        (∫_slow^shigh)(∫_tlow^thigh)  
            ElectField[x[s, t], y[s, t], z[s, t]] . normal[s, t] 
              dt ds
where normal[s, t] is as calculated above.

The flux is:

(∫_slow^shigh)
        (∫_tlow
          ^thigh) ElectField[x[s, t], y[s, t], z[s, t]] . normal[s, t] 
          \[DifferentialD]t \[DifferentialD]s

The flux measurement is 12  π flowing across the sphere from inside to outside.

-> To measure the flux through (x^2) + (y^2) + (z^2) = (r^2):
Use spherical coordinates for the sphere (x^2) + (y^2) + (z^2) = (r^2):

\(Clear[x, y, z, r, s, t, normal]; 
    
{x[s_, t_], y[s_, t_], z[s_, t_]} = {r\ Sin[s]\ Cos[t], 
        r\ Sin[s]\ Sin[t], r\ Cos[s]}; 

    \({{slow, shigh}, {tlow, thigh}} = {{0, π}, {0, 2\ π}}; \)\)\(normal[s_, t_] = 
      TrigExpand[
        D[{x[s, t], y[s, t], z[s, t]}, s]
            \[Cross]D[{x[s, t], y[s, t], z[s, t]}, t]]; 

    flux = ∫\_slow\%shigh
          \(∫\_tlow
            \%thigh ElectField[x[s, t], y[s, t], z[s, t]] . normal[s, t] 
            \[DifferentialD]t \[DifferentialD]s\)\)

This is the same as 12  π because the r terms cancel out.
This reveals that the flux across the sphere is 
       12  π 
no matter what the radius is.
Think of it this way:
The charge at {0, 0, 0} continually sends out barrages of "bullets" in all directions. That the flux across any sphere 
       (x^2) + (y^2) + (z^2) = (r^2)
is the same no mattter what r is means that at all times the same number of "bullets" passes through each of the spheres.

This also reflects the fact that this electric field has no sources or sinks other than the source at the singularity at {0, 0, 0}.  

All this is a direct consequence of the inverse square law.