VC.11 3D Flow Along
 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) Measuring flow along a 3D curve

B.1.a)

How do you measure the flow of a 3-dimensional vector field 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]} 
along a curve C given in parametric form
       P[t] = {x[t], y[t], z[t]} 
with a ≤ t ≤ b?

Just about the same way you do it in 2D.  
You measure the net flow of Field[x, y, z] along C by calculating the 3D path integral
       (∫_( C)^) 
        Field . unittan ds
        = 
        (∫_( C)^) m[x, y, z] dx +  
          n[x, y, z] dy +  
          p[x, y, z] dz
       =(∫_( C)^) (m[x[t], y[t], z[t]]x'[t]
		        + n[x[t], y[t], z[t]]y'[t]
            + p[x[t], y[t], z[t]]z'[t])dt
If this path integral calculates out as positive, then the net flow of Field[x, y, z] along C is in the direction specified by the parameterization.
If this path integral calculates out as negative, then the net flow of Field[x, y, z] along C is against the direction specified by the parameterization.

B.1.b)

Here is a 3D vector field:

\(Clear[Field, m, n, p, x, y, z]; 
m[x_, y_, z_] = y\^2\ z\^3; 

    \(n[x_, y_, z_] = 1\/3\ x\ y\ z\^3; \)\)p[x_, y_, z_] = (1/2) x (y^2) (z^2); Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}

Here's a curve C in three dimensions:

\(Clear[t]; 
{x[t_], y[t_], z[t_]} = {t, t + Sin[1.6\ t], t - 1}\)

This curve is given to start at:

a = 0.6; start = {x[a], y[a], z[a]}

And this curve is given to end at:

b = 2.8; end = {x[b], y[b], z[b]}

Here's a look at C:

\(Clear[P, t, Cplot, labels, setup]; 
P[t_] = {x[t], y[t], z[t]}; 

    \(Cplot = 
      ParametricPlot3D[Evaluate[P[t]], {t, a, b}, 
        DisplayFunction -> Identity]; \)\)\(labels = {Graphics3D[Text["start", start]], 
        Graphics3D[Text["end", end]]}; 


    setup = Show[Cplot, labels, ViewPoint -> CMView, 
        Axes -> Automatic, 
        AxesLabel -> {"x", "y", "z"}, Boxed -> False, 
        DisplayFunction -> $DisplayFunction]; \)

Measure the net flow of Field[x, y, z] along C, and interpret the result.

The measurement
       (∫_( C)^) 
        Field . unittan ds
        = 
        (∫_( C)^) m[x, y, z] dx +  
          n[x, y, z] dy +  
          p[x, y, z] dz
       =(∫_( C)^) (m[x[t], y[t], z[t]]x'[t]
		        + n[x[t], y[t], z[t]]y'[t]
            + p[x[t], y[t], z[t]]z'[t])dt
measures the flow of the vector field 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]}
along C:

NIntegrate[m[x[t], y[t], z[t]] x'[t]+n[x[t], y[t], z[t]] y'[t]+p[x[t], y[t], z[t]] z'[t],{t, a, b}]

Positive.
This means that the net flow of Field[x, y, z] along C is fairly strong in the direction of the parameterization of C.  
To see what this means, just plot C and some of its unit tangent vectors.

Clear[unittan];
unittan[t_]=TrigExpand[((P'[t])/(Sqrt[P'[t].P'[t]]))];
jump = ((b - a)/5);
Show[setup, 
        Table[Arrow[unittan[t], Tail -> P[t]], {t, a, b - jump, jump}]];

There you go with the flow.

B.2) The curl of a 3D vector field

Start with a cleared 3D vector field:

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

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

The curl, curlField[x, y, z], of Field[x, y, z] is given by:

\(Clear[curlField]; 

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

That's quite a pill to swallow, but the computer will swallow it for you.
You'll want to use it to help finger a 3D vector field.

B.2.a) Using curlField[x, y, z] to help finger a 3D vector field.

Stick the tail of a unit vector V at a point {a, b, c}.
Push your finger onto the point of V so that V spikes through the center of the bones of your finger, and the tip of your finger is at {a, b, c}. 
How do you use 
       curlField[a, b, c] 
to determine whether the tip of your finger feels a net clockwise or counterclockwise swirl resulting from the flow of the given 3D vector field
       Field[x, y, z]?

Here's how it works:
-> If curlField[a, b, c] . V > 0,
then you feel Field[x, y, z] swirling in the counterclockwise way at the tip of your finger.
-> If curlField[a, b, c] . V < 0,
then you feel Field[x, y, z] swirling in the clockwise way at the tip of your finger.
-> If curlField[a, b, c] . V = 0, 
then you feel no net swirl at all.
Try it out for a given vector field, a given unit vector V, and a given point {a, b, c}:

Clear[Field, curlField, unitvector, s, t, x, y, z, m, n, p]; unitvector[s_, t_] := N[{Sin[s] Cos[t], Sin[s] Sin[t], Cos[s]}]m[x_, y_, z_] = x - 2 z; n[x_, y_, z_] = y - 2 (x^2); p[x_, y_, z_] = z - 2 y; Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}; \(curlField[x_, y_, z_] = {D[p[x, y, z], y] - D[n[x, y, z], z], 
        D[m[x, y, z], z] - D[p[x, y, z], x], 
        D[n[x, y, z], x] - D[m[x, y, z], y]}; 

    curlField[0, 2, 1] . unitvector[π\/3, π\/4]\)

Negative.
This tells you that when you push your finger onto 
       V = unitvector[(π/3), (π/4)] 
so that V spikes through the center of your finger and the tip of your finger is at {0, 2, 1}, then your finger feels a net CLOCKWISE swirl.
Confirm this by plotting what's happening near {0, 2, 1} using  the plotting option ViewPoint -> V to look from the tip of V to the tail of V at the point {0, 2, 1}:

point = {0, 2, 1}; V = unitvector[(π/3), (π/4)]; h = 0.3; scalefactor = 0.08; \(fieldplot = 
      Table[Arrow[Field@@\((point + {x, y, z})\), 
          Tail -> N[point + {x, y, z}], VectorColor -> Blue, 
          ScaleFactor -> scalefactor], {x, \(-h\), h, 2\ h}, {y, \(-h\), 
          h, 2\ h}, {z, \(-h\), h, 2\ h}]; 

    pointplot = Graphics3D[{PointSize[0.07], Point[point]}]; 


    swirlplot = 
      Show[pointplot, fieldplot, PlotRange -> All, Boxed -> False, 
        ViewPoint -> V]; \)

Put the tip of your finger at the plotted point on the screen and feel that net clockwise flow!
Just as the calculation predicted!
Try it for a different point and a different unit vector:

curlField[0, 0, 0] . unitvector[((3 π)/2), (π/6)]

Positive. 
This tells you that when you push your finger onto 
       V = unitvector[((3 π)/2), (π/6)] 
so that V spikes though the center of your finger and the tip of your finger is at {0, 0, 0}, then your finger feels a net COUNTERCLOCKWISE swirl.
Confirm with a plot:

point = {0, 0, 0}; V = unitvector[((3 π)/2), (π/6)]; h = 0.2; scalefactor = 0.4; \(fieldplot = 
      Table[Arrow[Field@@\((point + {x, y, z})\), 
          Tail -> N[point + {x, y, z}], VectorColor -> Blue, 
          ScaleFactor -> scalefactor], {x, \(-h\), h, 2\ h}, {y, \(-h\), 
          h, 2\ h}, {z, \(-h\), h, 2\ h}]; 

    pointplot = Graphics3D[{PointSize[0.07], Point[point]}]; 


    swirlplot = 
      Show[pointplot, fieldplot, PlotRange -> All, Boxed -> False, 
        ViewPoint -> V]; \)

Put the tip of your finger at the plotted point on the screen and feel that net clockwise flow! 
Math happens again.

B.3) Stokes's formula for using the curl to measure the swirl of a vector field in 3D

In 2D, everyone likes the form of the Gauss-Green formula that tells you that if 
       Field[x, y] = {m[x, y], n[x, y]} 
and R is a region in two dimensions with boundary curve C, then 
       (∫_^)(∫_( R)) rotFielddA
        = 
        (∮_( C)) 
          Field . unittan ds
        = 
        (∮_( C)) m[x, y] dx  +  
          n[x, y] dy
        = net flow of Field[x, y] along C.
Stokes's formula is the 3D analogue of this formula.  Stokes's 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 R is a surface in three dimensions with boundary curve C, then
       (∫_^)(∫_( R)^) curlField.topunitnormaldA
        = 
        (∮_( C)) 
          Field . unittan ds
        = 
        (∮_( C)) m[x, y, z] dx  +  
          n[x, y, z] dy +  
          p[x, y, z] dz
        = net flow of Field[x, y, z] along C.
Before you can make complete sense of this formula, you have to deal with a few objections.

B.3.a.i) Two objections.

-> Objection 1:
The integral 
       (∮_( C)) Field . unittan ds
insists on a counterclockwise parameterization for the boundary curve C.  But counterclockwise parameterization makes no sense in three dimensions, so Stokes's formula  
       (∫_^)(∫_( R)^) curlField.topunitnormaldA = (∮_( C)) 
                Field . unittan ds
makes no sense.
-> Objection 2:
Look at this piece of the yz-plane:

yzplanepiece = 
      Graphics3D[Polygon[{{0, 0, 0}, {0, 2, 0}, {0, 2, 1}, {0, 0, 1}}]]; 
    Show[yzplanepiece, Axes -> True, 
      AxesLabel -> {"x", "y", "z"}, ViewPoint -> CMView, 
      Boxed -> False, PlotRange -> All, 
      DisplayFunction -> $DisplayFunction];

Because this surface sits in the yz-plane, you can't say whether the front side or the back side is the top side.
Calculating
       (∫_^)(∫_( R)^) curlField.topunitnormaldA
depends on knowing which side of the surface is the top side. So again Stokes's formula  
       (∫_^)(∫_( R)^) curlField.topunitnormaldA = (∮_( C)) 
                Field . unittan ds
makes no sense.
How do you get around these legitimate objections?

Taken separately, each of these two objections ruins any possibility of making sense out of Stokes's formula, but if you deal with both of them simultaneously, then you can get a full understanding of what Stokes's formula says.
First, you pick one side of the surface at your own pleasure and call it the top side.

Once you pick the top side of R, then you agree that you are going in the counterclockwise direction on C if, as you walk around C, R is on your left.  With these two agreements carved in stone, Stokes's formula       
       (∫_^)(∫_( R)^) curlField.topunitnormaldA = (∮_( C)) 
                Field . unittan ds
       ( = 
          (∫_( C)^) m[x, y, z] dx +  
            n[x, y, z] dy +  p[x, y, z] dz)
        
for a given 3D vector field 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]} 
and a given surface R with boundary curve C makes perfect sense provided you agree to parameterize C in the counterclockwise direction dictated by your choice of the top side of the surface according to the agreement above.

B.3.a.ii) The designated top side.

Here is a surface R which is a cap of an ellipsoid:

\(Clear[x, y, z, s, t]; 
    
{x[s_, t_], y[s_, t_], z[s_, t_]} = {5\ Sin[s]\ Cos[t], 
        4\ Sin[s]\ Sin[t], 6\ Cos[s]}; 
\({slow, shigh} = {0, π\/3}; 
    \)\)surfaceplot = 
      ParametricPlot3D[{x[s, t], y[s, t], z[s, t]}, {s, slow, shigh}, {t, 0, 
          2 π}, Boxed -> False, PlotRange -> All, 
        ViewPoint -> CMView, 
        AxesLabel -> {"x", "y", "z"}];

A parameterization of the boundary curve C is:

\(Clear[P]; 
P[t_] = {x[shigh, t], y[shigh, t], z[shigh, t]}\)

Here are the surface R, its boundary curve C, and a few tangent vectors reflecting the direction of this parameterization of C:

Cplot = 
      ParametricPlot3D[Evaluate[P[t]], {t, 0, 2 π}, 
        DisplayFunction -> Identity]; tangents=Table[Arrow[P'[t],Tail -> P[t]],{t, 0, 2 π, (π/3)}];Show[surfaceplot, Cplot, tangents, Boxed -> False, 
      PlotRange -> All, AxesLabel -> {"x", "y", "z"}, 
      ViewPoint -> CMView, DisplayFunction -> $DisplayFunction];

If you want to call this parameterization of boundary curve C counterclockwise in accordance with the agreement above, then which side of the cap must you designate as the top side?

Walk around the boundary in the direction indicated by the tangents.  As you walk, your left foot hits on the visible side (the high side). This is your designated top side.

B.3.a.iii) A third objection:

Take a look at this sorry excuse for a surface:

Clear[P, x, y, z, t, unittan, mainunitnormal, binormal, moebius];
P[t_] = 2 {-Cos[t], Sin[t], 0};
unittan[t_]=TrigExpand[((P'[t])/(Sqrt[P'[t].P'[t]]))];
mainunitnormal[t_]=unittan'[t]/(√TrigExpand[unittan'[t].unittan'[t]]);
binormal[t_] = TrigExpand[unittan[t]\[Cross]mainunitnormal[t]];
moebius[s_, t_] = 
          P[t] + s 
              (Sin[(t/2)] mainunitnormal[t] - Cos[(t/2)] binormal[t]);
ParametricPlot3D[
          Evaluate[moebius[s, t]], {s, -0.5, 0.5}, {t, -π, 
            π}, PlotPoints -> {2, 40}, ViewPoint -> CMView, 
          PlotRange -> All, Boxed -> False, Axes -> False];

This surface, which folks like to call a "Moebius strip," is unfortunate enough to have only one side.  How the hell are you going to designate a top side for a surface that has only one side?

You can't. 
It's impossible to choose a top side.
The upshot:
It's impossible to use Stokes's formula for a surface that has only one side.

B.3.a.iv)

What else do you have to worry about?

Stokes's formula works only for surfaces that are not the complete boundary of a solid region. 
For example, Stokes's formula is OK if R is a cap of a sphere, but Stokes's formula does not work if R is the whole sphere.  This is the case simply because a whole sphere does not have a boundary curve, but any cap of a sphere does have a boundary curve.

B.3.b.i) Stokes's formula helps to explain why the finger test works.

Take a 3D vector field Field[x, y, z].
Stick the tail of a unit vector V at the point {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}.
Push your finger onto V so that V spikes through the center of your finger and the tip of your finger is at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}.
Why does the tip of your finger feel a counterclockwise swirl if 
       curlField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] . V > 0?
Why does the tip of your finger feel a clockwise swirl if 
       curlField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] . V < 0?

Put a plane with normal vector V through the point 
       {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}. 
Put a small circle C with center at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])} in this plane. Call R the part of this plane that is within the small circle. Designate the top side by taking the direction V points.  This gives you
       V = topunitnormal
for this little surface R. 
Stokes's formula tells you:
       net flow of Field[x, y, z] along C 
        = 
        (∫_( C)^) m[x, y, z] dx +  
          n[x, y, z] dy +  
          p[x, y, z] dz
        = 
        (∫_^)(∫_( R)^) 
            curlField . topunitnormal dA
        = 
        (∫_^)(∫_( R)^) 
            curlField . V dA.
Here's the kicker:  
If 
       curlField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] . V > 0, 
then  
       curlField[x, y, z] . V > 0 
for all {x, y, z}'s close to {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}.  So if C is small enough so that 
       curlField[x, y, z] . V > 0
at all {x, y, z}'s on R and inside C, then you see that
       net flow of Field[x, y, z] along C 
       =(∫_^)(∫_( R)^) curlField.VdA>0.
So the net flow of Field[x, y, z] along such a small circle C is counterclockwise.
This is why the tip of your finger feels a counterclockwise swirl when you stick the tail of V at the point {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}, and push your finger onto V so that V spikes through the center of your finger and the tip of your finger is at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}.
 
Similarly, if curlField[(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])] . V < 0 and you stick the tail V at the point {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}, and push your finger onto V so that V spikes through the center of your finger and the tip of your finger is at {(Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])}, then your fingertip feels a clockwise swirl.

B.3.b.ii)

Given a 3D vector field Field[x, y, z], when you stick the tail of a unit vector V at the point {x, y, z} and push your finger onto V so that V spikes through the center of your finger and the tip of your finger is at {x, y, z}, then what unit vector V do you want to use so that your fingertip feels the greatest possible counterclockwise swirl?
What does the direction of curlField[x, y, z] tell you?
What does the length of curlField[x, y, z] tell you?

The swirl your fingertip feels is measured by
       curlField[x, y, z] . V 
        = 
        ‖curlField[x, y, z]
              ‖ 
          ‖V‖ 
          Cos[angle between]
        = 
        ‖curlField[x, y, z]
              ‖ Cos[angle between]
because V is a unit vector.
This is biggest when 
       Cos[angle between] = 1 = Cos[0].
So you feel the biggest counterclockwise swirl when you take the unit vector V to be in the same direction as curlField[x, y, z].
In other words, to feel the most counterclockwise swirl, you take 
       V=curlField[x, y, z]
                /(Sqrt[curlField[x, y, z] . curlField[x, y, z]])
The upshots: 
-> The direction of curlField[x, y, z] gives the direction of the axis of the actual swirl of Field[x, y, z] at {x, y, z}.
-> The length of curlField[x, y, z] measures how vigorous the swirl at {x, y, z} is.
Check this out:

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

    m[x_, y_, z_] = x + 2\ z\^2; 
\(n[x_, y_, z_] = y - 2\ x\^2; \)\)p[x_, y_, z_] = z + 4 (y^2); Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}; curlField[x_, y_, z_] = {D[p[x, y, z], y] - D[n[x, y, z], z], 
        D[m[x, y, z], z] - D[p[x, y, z], x], 
        D[n[x, y, z], x] - D[m[x, y, z], y]}

To feel the greatest possible counterclockwise swirl of this vector field at {1, 1, 1}, you should use:

point = {1, 1, 1}; V = 
      curlField[1, 1, 1]/(Sqrt[curlField[1, 1, 1] . curlField[1, 1, 1]])

Take a look down the vector V from its tip to its tail at {1, 1, 1}:

h = 0.3; scalefactor = 0.1; \(fieldplot = 
      Table[Arrow[N[Field@@\((point + {x, y, z})\)], 
          Tail -> N[point + {x, y, z}], VectorColor -> Blue, 
          ScaleFactor -> scalefactor], {x, \(-h\), h, 2\ h}, {y, \(-h\), 
          h, 2\ h}, {z, \(-h\), h, 2\ h}]; 

    pointplot = Graphics3D[{PointSize[0.07], Point[point]}]; 


    fingerplot = 
      Show[pointplot, fieldplot, Boxed -> False, ViewPoint -> V]; 
    \)

Powerfully counterclockwise just as you expected.
Math really happens.

B.3.b.iii)

Where did the name "curl" come from?

Take a look at the end of the answer to part ii) above.

B.3.c.i)

Take a surface R with boundary curve C.  Designate a top side of R.  The path integral  
       (∫_( C)^) 
        Field . unittan ds
        = 
        (∫_( C)^) m[x, y, z] dx +  
          n[x, y, z] dy +  
          p[x, y, z] dz
measures the net flow of a given 3D vector field 
       Field[x, y, z] = {m[x, y, z], n[x, y, z], p[x, y, z]} 
along C. 
If you know that
       curlField[x, y, z] . topunitnormal[x, y, z] > 0 
at all points {x, y, z} on the surface R, then how do you also know that the net flow of Field[x, y, z] along the boundary curve C is counterclockwise?

You calculate the flow of Field[x, y, z] along C by calculating 
       (∫_( C)^) m[x, y, z] dx +  
        n[x, y, z] dy +  p[x, y, z] dz.
Stokes's formula tells you that this is the same as calculating
       (∫_^)(∫_( R)^) 
          curlField . topunitnormal dA
Here's the kicker:  
If 
       curlField[x, y, z] . topunitnormal[x, y, z] > 0 
at all points of R, then the 
       flow of Field[x, y, z] along C 
        = 
        (∫_( C)^) m[x, y, z] dx +  
          n[x, y, z] dy +  
          p[x, y, z] dz
        = 
        (∫_^)(∫_( R)^) 
              curlField . topunitnormal dA > 0,
because curlField[x, y, z] . topunitnormal[x, y, z] > 0 for all {x, y, z}'s on R.
So, if you know that 
       curlField[x, y, z] . topunitnormal[x, y, z] > 0
for all {x, y, z}'s on R, then you know that the flow - along - C measurement is positive.  
This means that the net flow of Field[x, y, z] along the boundary curve C of R is counterclockwise.

B.3.c.ii)

Take a surface R with boundary curve C. If you know that
       curlField[x, y, z] = 0 
at all points {x, y, z} on the surface R, then how do you also know that the net flow of Field[x, y, z] along C is 0?

You calculate the net flow of Field[x, y, z] along C by calculating 
       (∫_( C)^) m[x, y, z] dx +  
        n[x, y, z] dy +  p[x, y, z] dz.
Stokes's formula tells you that this is the same as calculating
       (∫_^)(∫_( R)^) 
          curlField . topunitnormal dA.
Here's the kicker:  
If 
       curlField[x, y, z] = 0 
at all points of R, then the 
       net flow of Field[x, y, z] along C 
        = 
        (∫_( C)^) m[x, y, z] dx +  
          n[x, y, z] dy +  
          p[x, y, z] dz
        = 
        (∫_^)(∫_( R)^) 
            curlField . topunitnormal dA
        = 
        (∫_^)(∫_( R)^) 
            0 dA
        = 0, 
because curlField[x, y, z] = 0 for all {x, y, z}'s on R.
So 
       curlField[x, y, z] = 0 
for all {x, y, z}'s on R guarantees that the flow - along - C measurement is 0.  This means that the net flow of Field[x, y, z] along the boundary curve C of R is 0.

B.3.c.iii) Irrotational vector fields.

Why do the fancy folks call a vector field Field[x, y, z] with 
       curlField[x, y, z] = 0 
at all points {x, y, z} irrotational?

As you saw above if 
       curlField[x, y, z] = 0 
for all points {x, y, z}, then the flow of Field[x, y, z] along any closed curve is 0. 
This rules out any vortex or other swirl.

B.4) An attempt to explain Stokes's formula as an outgrowth of the 2D Gauss-Green formula

B.4.a)

Explain the reasoning behind Stokes's formula.

Ultimately Stoke's formula is a consequence of the Gauss-Green formula in two dimensions which, in turn, is a consequence of the fundamental formula of calculus in one variable. 
The fundamental formula of calculus strikes again!
You can explain Stokes's formula as a rather easy consequence of the Gauss-Green formula, but there is a lot of bookkeeping (bean counting).
The first step is to reduce the explanation of the overall formula to an explanation for a very simple surface.

Step One:  Reduction to parallelograms.

Look at a picture and think of it as a view of a surface R in three dimensions.

surface = 
      Plot[{(Sqrt[1 - (x^4)]), -(Sqrt[1 - (x^4)])}, {x, -1, 1}, 
        AspectRatio -> Automatic, Axes -> None];

Split the surface into two pieces and label:

splitter = 
      Graphics[{RGBColor[1, 0, 0], Line[{{0, 1}, {0, -1}}]}]; labels = 
      Graphics[{Text["Part1", {-(1/2), 0}], Text["Part2", {(1/2), 0}], 
          Text["edge1", {1, 0}], Text["edge2", { .15, - .2}], 
          Text["edge3", {- .15,  .2}], Text["edge4", {-1, 0}]}]; 
    \(arrows = {Vector[{ .9, \(-\(1\/3\)\)}, { .9, 1\/3}], 
        Vector[{\(- .9\), 1\/3}, {\(- .9\), \(-\(1\/3\)\)}], 
        Vector[{ .1, 1\/3}, { .1, \(-\(1\/3\)\)}], 
        Vector[{\(- .1\), \(-\(1\/3\)\)}, {\(- .1\), 1\/3}]}; 

    \(Show[{surface, splitter, labels, arrows}, 
      PlotRange -> {{\(-1.2\), 1.2}, {\(-1.1\), 1.1}}]; \)\)

Now instead of one surface, you've got two surfaces with their own boundaries.  To integrate a function around the boundary of Part2, you integrate out around edge1, and then down along edge2.  To integrate around the boundary of Part1, you integrate up along edge3 and out around edge4.  
Look at what happens:
       (∫_( edge1)^) F . unittan ds +  
        (∫_( edge2)^) F . unittan ds
               
        + (∫_( edge3)^) 
            F . unittan ds +  
        (∫_( edge4)^) 
          F . unittan ds
        = 
        (∫_( edge1)^) F . unittan ds +  
          (∫_( edge4)^) 
            F . unittan ds,
because edge2 and edge3 are the same curves going in opposite directions.   
What does it say?  
It says that adding the integrals around the boundaries of the surfaces Part1 and Part2 gives the same result as integrating around the boundary of the whole region R before you split it.  
And what does that say to you? 
It tells you that to explain Stokes's formula   
       (∫_^)(∫_( R)^) 
            curlField . topunitnormal dA = 
        (∮_( C)^) 
          Field . unittan ds,
you can cut the original surface R up into small pieces, and if you can explain the formula for the small pieces and their boundary curves, then you can add the formulas for the small pieces to arrive at a full explanation of the formula for the whole original surface R.
Here is another view of R but this time R is cut up into little pieces:

Plot3D[
      (Sqrt[1 - (x^2) - (y^2)]), {x, -(1/2), (1/2)}, {y, -(1/2), 
        (1/2)}, Boxed -> False, Axes -> None];

So explaining the formula for the whole surface is the same as justifying it for each of the little parallelogram-like parts you see above and adding the results together.  
Now things are a lot simpler because each of the smaller parts is very close to being a parallelogram.  So if you can explain why Stokes's formula holds for a parallelogram, you'll have explained Stokes's formula for R.

Step Two: Explaining Stokes's formula for a parallelogram.

You're not out of the woods yet.  This is the part that is  complicated ... not hard, just complicated.  It's really a job for those bean counting dweebs who like to do a lot of bookkeeping.  Recall what needs to be explained: 
       (∫_^)(∫_( R)^) 
            curlField . topunitnormal dA = 
        (∮_( C)^) 
          Field . unittan ds,
where the surface R is a parallelogram, and its boundary C is the edge of the parallelogram.  Here's a picture with normal vector included:

a = (1/3) {1, 1, 1}; aa = (1/3) {1, 1, 0}; b = {-(1/2), 0, (1/4)}; bb = {-(1/2), 0, 0}; c = {0, -(1/3), (1/8)}; cc = {0, -(1/3), 0}; para = 
      Graphics3D[{RGBColor[1, 0.5, 0], 
          Polygon[{a, a + b, a + b + c, a + c, a}]}]; base = 
      Graphics3D[{RGBColor[0.8, 0.8, 0.8], 
          Polygon[{aa, aa + bb, aa + bb + cc, aa + cc, aa}]}]; \(normal = Vector[a + \(b + c\)\/3, a + \(b + c\)\/3 + 2\ b\[Cross]c]; 

    lines = Graphics3D[{Line[{a, aa}], Line[{a + b, aa + bb}], 
          Line[{a + c, aa + cc}], Line[{a + b + c, aa + bb + cc}]}]; 


    Show[lines, normal, para, base, ViewPoint -> CMView, 
      Ticks -> None, Axes -> Automatic, 
      AxesLabel -> {"x", "y", "z"}, Boxed -> False, 
      PlotRange -> All]; \)

The upper parallelogram is the region R you are working with, and the lower parallelogram is the shadow of the upper parallelogram on the xy-plane.  When you parameterize to integrate around the boundary C on the top, you will be directly over a counterclockwise parameterization of the boundary of the lower parallelogram in the xy-plane.  Ultimately you are going to use the two-dimensional Gauss-Green formula applied to the lower parallelogram which will be called (Subscript[R, uv]) when you parameterize with x = u and y = v.  The boundary of (Subscript[R, uv]) will be called (Subscript[C, uv]).
Now set up some calculations:

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

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

The upper parallelogram with the normal vector pointing up out of it lies in a plane.  Write the equation for this plane as 
       z = c +  a x +  b y, 
or parametrically by:

\(Clear[u, v, a, b, c]; 
    
{x[u_, v_], y[u_, v_], z[u_, v_]} = {u, v, c + a\ x[u, v] + b\ y[u, v]}
        \)

As {u, v} varies through (Subscript[R, uv]), {x[u, v], y[u, v], z[u, v]} varies through R.
Look at the normal vector:

\(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]\)

Because this normal has a positive third slot, it points up, just as in the picture above.

\(Clear[topunitnormal]; 

    topunitnormal[u_, v_] = 
      normal[u, v]\/\@\(normal[u, v] . normal[u, v]\)\)

This means that when you calculate the path integral you will have to parameterize the boundary curve of the upper parallelogram so that as you walk in the direction of the parameterization, your left foot touches the high side of the surface R.
Now see what 	
       (∫_^)(∫_( R)^) 
          curlField . topunitnormal dA
        = 
        (∫_^)(∫_( (Subscript[R, uv]))^) 
            curlField[x[u, v], y[u, v], z[u, v]] . topunitnormal[u, v] 
              (Subscript[SA, xyz])[u, v] du dv   looks like:

\(Clear[SAxyz]; 
SAxyz[u_, v_] = \@\(normal[u, v] . normal[u, v]\)\)

\(Clear[curlField]; 

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

The integrand is:

integrand = 
      Simplify[curlField[x[u, v], y[u, v], z[u, v]] . topunitnormal[u, v] 
          SAxyz[u, v]]

Here's the bookkeeping trick:  
Once you have done this, you can reduce the whole explanation to the two dimensional Gauss-Green formula.  Start by making a couple of off-the-wall but delicious choices for new 2-dimensional functions.

\(Clear[mm, nn, FF, rotFF]; 

    mm[u_, v_] = 
      m[x[u, u], y[u, v], z[u, v]] + a\ p[x[u, v], y[u, v], z[u, v]]; 

    \(nn[x_, y_] = 
      n[x[u, v], y[u, v], z[u, v]] + b\ p[x[u, v], y[u, v], z[u, v]]; \)\)FF[u_, v_] = {mm[u, v], nn[u, v]}

Look at rotFF[u, v]:

rotFF[u_, v_] = Expand[D[nn[u, v], u] - D[mm[u, v], v]]

rotFF[u, v] looks like the quantity you got for the integrand up above.  
Is it?

integrand == rotFF[u, v]

You bet it is!    
This tells you that
       (∫_^)(∫_( R)^) 
          curlField . topunitnormal dA
        = 
        (∫_^)(∫_( (Subscript[R, uv]))^) 
            rotFF[u, v] du dv
and the two dimensional Gauss-Green formula says that
       (∫_^)(∫_( (Subscript[R, uv]))^) 
          rotFF[u, v] du dv
        = 
        (∮_( (Subscript[C, uv]))^) mm[u, v] du + 
           nn[u, v] dv
where (Subscript[C, uv]) is the boundary of (Subscript[R, uv]) .
Now 
       (∮_( (Subscript[C, uv]))^) mm[u, v] du + 
         nn[u, v] dv        = 
        (∮_( (Subscript[C, uv]))^) 
            (m[x[u, v], y[u, v], z[u, v]] +  
                  a p[x[u, v], y[u, v], z[u, v]]) du  +  
          (n[x[u, v], y[u, v], z[u, v]] +  b p[x[u, v], y[u, v], z[u, v]])
               dv    

 = 
        (∮_( (Subscript[C, uv]))^) 
            m[x[u, v], y[u, v], z[u, v]] du +  
          n[x[u, v], y[u, v], z[u, v]] dv +  
          p[x[u, v], y[u, v], z[u, v]] 
            (a du +  b dv) . 
Now look at:

z[u, v]

Read off 
       dz = 
        a du +  b dv, 
and plug into the above, remembering that x = u and y = v, to realize that when you parameterize and calculate
       (∮_( C)^) 
        Field . unittan ds
        = 
        (∮_( C )^) m[x, y, z] 
              dx  +  n[x, y, z] dy  +  
          p[x, y, z] dz, 
then you get 
       (∮_( (Subscript[C, uv]))^) mm[u, v] du + 
         nn[u, v] dv
        = 
        (∫_^)(∫_( (Subscript[R, uv]))^) 
            rotFF[u, v] du dv
        = 
        (∫_^)(∫_( R)^) 
            curlField . topunitnormal dA.
And now the formula
       (∮_( C)^) 
          Field . unittan ds = (∫_^)(∫_( R)^) 
          curlField . topunitnormal dA
has been explained.
Review the explanation to see why Stokes's formula is just a dressed-up version of the basic two dimensional Gauss-Green formula.