How can I use ContourPlot to creat streamline: proper arrow orientation and correct domain?

Question

I get into trouble when I use ContourPlot to draw streamlines in a free surface flow domain.

This is my data.The data are a little bit big, please be patient. Here is my code:

Get["C:\\Users\\Administrator\\Desktop\\...\\data.mx"]

Firstly, I plot the free surface, which is a function of x.

profile = Plot[mdfun[x, 6], {x, 0, 2*Sqrt[2] \[Pi]}, 
PlotRange -> {{0, 2*Sqrt[2] \[Pi]}, {0, 2}},
ImageSize -> 400,Frame -> True, Axes -> False,
PlotStyle -> {Directive[Black, Thick]},
FrameTicks -> {{{0, 0.2, 0.4, 0.6, 0.8, {1, "1.0"}, 1.2, 1.4, 1.6, 1.8}, None},
{{0, {1, ""}, 2, {3, ""}, 4, {5, ""}, 6, {7, ""}, 8}, None}}]

hx = Derivative[1, 0][mdfun]; (*get the first derivative with respect to x*)
hxxx = Derivative[3, 0][mdfun];  (*get the third derivative with respect to x*)
scaledA[d_, B_, k_, e_, hx_, h_, hxxx_] := (-3*d*(1 + B*k) *hx )/((k + h + B*k*h)^3*e) - 3*1/e*hxxx - 3*1/e*hx;
A = scaledA[1, 1, 0.1, 0.05, hx[x, 6], mdfun[x, 6], hxxx[x, 6]]; (*calculate a term of stream function*)
scaled\[Phi][A_, z_, h_, B_, k_, hx_, e_, m_] := A*(z^3/6 - (h*z^2)/2) + m ((1 + B k)*k*hx*z^2)/((k + h + B k h)^2*e); (* stream function *)

sl = ContourPlot[scaled\[Phi][A, y, mdfun[x, 6], 1, 0.1, hx[x, 6], 0.05, 1] == Range[0.1, 1.3, 0.2],{x, 0, 2*Sqrt[2] \[Pi]}, {y, 0, mdfun[x, 6]}, ImageSize -> 400,ContourStyle -> {
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin]
}, Frame -> False (*,PerformanceGoal\[Rule]"Speed"*)] (*plot seven streamlines *)

NOTE: ContourPlot takes about 60s on a common PC, you can uncomment PerformanceGoal\[Rule]"Speed" to speed up. Please be patient...The reason why I plot 7 streamlines is to show the different directions of the arrows drawn

Then, I got the streamlines

If I use RegionFunction to explicitly specify the domain drawn.

sl2 = ContourPlot[scaled\[Phi][A, y, mdfun[x, 6], 1, 0.1, hx[x, 6], 0.05, 1] == Range[0.1, 1.3, 0.2],
{x, 0, 2*Sqrt[2] \[Pi]}, {y, 0, mdfun[x, 6]}, ImageSize -> 400, 
RegionFunction -> Function[{x, y, z}, 0 <= y <= mdfun[x, 6]],
ContourStyle -> {
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin],
Directive[Black, Thin]
}, Frame -> False(*,PerformanceGoal\[Rule]"Speed"*)]

I still get the same plot as above.

Then I add arrow for the streamlines

asl = sl /. Line[x_] :> {Arrowheads[{{0.02, 0.2}, {0.02, 0.4}, {0.02, 0.6}, {0.02, 0.85}}], Arrow[x]}

And I got the following plot

Last, I combine the profile and streamlines.

Show[profile, asl]

Here I have four questions:

(1) How can I specify the upper limit in a free surface problem? Here I employed {y, 0, mdfun[x, 6]};

(2) Is there any method to find all the representative or distinctive streamlines instead of assign different values to ContourPlot every time;

(3) Frankly speaking, I am not sure whether my streamline is correct or has the completed range. Because I drew them by restricting the y-range by {y, 0, mdfun[x, 6]} in the first code, in addition, by RegionFunction -> Function[{x, y, z}, 0 <= y <= mdfun[x, 6]] in the second code sl2. I really do not understand the difference between them in my case.

(4)This question has been solved by @ bbgodfrey. As you can see, the orientation of the arrows of the streamlines is different from each other. I have tried to change the sign of 0.02 in Arrowheads, but I still can not get the same orientation anyway, say towards the right. Can any one tell me how can I manipulate the arrow individually for each streamline?

Thank you for your patience!

Answer 1

Although InputForm for these plots is enormous, due to the InterpolatingFunctions used,

Cases[asl, GraphicsComplex[x__]]

allows the data associated with Arrow and Arrowheads to be isolated. Doing so shows that the data for most contours runs right to left, but for a few from left to right.

A work-around is as follows. Determine the Position of the Line Commands in sl.

Position[sl, Line[x_], Infinity]
(*{{1, 2, 1, 3, 1, 2}, {1, 2, 1, 3, 1, 3}, {1, 2, 1, 3, 1, 4},
{1, 2, 1, 3, 1, 5}, {1, 2, 1, 3, 1, 6}, {1, 2, 1, 3, 1, 7}, {1, 2, 1, 3, 1, 8}}*)

The first entry is associated with the bottom curve, the next entry with the second curve, etc. Use this information to reverse the order of points in the desired lines. Here, I assume that arrow heads in lines 3 and 6 are to be reversed. So, insert

Part[sl, 1, 2, 1, 3, 1, 4, 1] = Reverse[sl[[1, 2, 1, 3, 1, 4, 1]]];
Part[sl, 1, 2, 1, 3, 1, 7, 1] = Reverse[sl[[1, 2, 1, 3, 1, 7, 1]]];

just before the line creating asl to produce

as desired. Inelegant but effective.

Addendum

The process described above can be automated using a command structurally similar to that used by Michael E2. Instead of reversing the data of curves on a case by case basis, replace the asl command by

asl = sl /. gc : GraphicsComplex[p_, __] :> (gc /. 
  Line[x_] :> {Arrowheads[{{0.02, 0.2}, {0.02, 0.4}, {0.02, 0.6}, {0.02, 0.85}}], 
  Arrow[If[p[[First@x, 1]] > p[[Last@x, 1]], x, Reverse[x]]]})

Answer 2

Along the lines of Rahul's comment:

Two possible methods:

1. If the function f is differentiable, the compare the directions of the Grad of f with the (2D) Cross of the direction of the line.

2. Use the Cross of the direction of the line to determine which side of the contour has the higher value of f by comparing the value of f on contour with a point on the side (toward which Cross points).

Toy example (I worked this out before seeing Rahul's comments):

f[x_, y_] := y - x^4/100 + x^2/10;
sl = ContourPlot[
  f[x, y] == Range[-1, 0.4, 0.2], {x, -4, 4}, {y, -2, 1}, 
  ImageSize -> 400, ContourStyle -> Directive[Black(*,Thin*)], 
  Frame -> False (*,
  PerformanceGoal\[Rule]"Speed"*)] (*plot seven streamlines*)

asl = sl /. 
  Line[x_] :> {Arrowheads[{{0.02, 0.2}, {0.02, 0.4}, {0.02, 0.6}, {0.02, 0.85}}], Arrow[x]}

Method 1:

asl = sl /. 
  gc : GraphicsComplex[p_, __] :> (gc /. 
     Line[ll_] :> {Arrowheads[{{0.02, 0.2}, {0.02, 0.4}, {0.02, 0.6}, {0.02, 0.85}}],
       Arrow[If[
         Dot[
           Grad[f[x, y], {x, y}] /. Thread[{x, y} -> p[[First@ll]]],
           Cross[First@Differences[p[[ll[[;; 2]]]]]]
           ] > 0,
         Reverse[ll], ll]]})

Method 2:

scale = 0.1; (* controls how far to step to the side *)
asl = sl /. 
  gc : GraphicsComplex[p_, __] :> (gc /. 
     Line[ll_] :> {Arrowheads[{{0.02, 0.2}, {0.02, 0.4}, {0.02, 0.6}, {0.02, 0.85}}],
       Arrow[If[
         Less @@ (f[x, y] /. {Thread[{x, y} -> p[[ll[[2]]]]],
             Thread[{x, y} -> 
               p[[First@ll]] + scale Cross[First@Differences[p[[ll[[;; 2]]]]]]
              ]}),
         Reverse[ll], ll]]})

Output:

The orientation can be reverse by reversing the comparison in the If statement.