9
$\begingroup$

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

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

enter image description here

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

enter image description here

Last, I combine the profile and streamlines.

Show[profile, asl]

enter image description here

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!

$\endgroup$
  • $\begingroup$ If all you're trying to do is to draw arrows in the contour lines, why did you post your full problem? Just make up a fake and simple CountourPlot[] and re-ask without all the burden $\endgroup$ – Dr. belisarius Dec 19 '14 at 15:48
  • $\begingroup$ Here's a nice fake example for y'all: sl = ContourPlot[Cos[x] + 1/2 y^2, {x, -2 π, 2 π}, {y, -π, π}, Contours -> 5, AspectRatio -> Automatic, ContourStyle -> Directive[Black, Opacity[1]], ContourShading -> None] $\endgroup$ – Rahul Dec 19 '14 at 16:45
  • 2
    $\begingroup$ Anyway, I do think that this is a misfeature of Mathematica. Contour lines of a scalar function can be assigned a consistent orientation based on which side of the line has the larger function value, and it would be nice if Mathematica automatically respected that orientation. In the meantime, I hope someone posts a nice workaround. $\endgroup$ – Rahul Dec 19 '14 at 16:48
  • $\begingroup$ @ belisarius. Please see the update, the questions are related to my full problem. Thanks :) $\endgroup$ – Enter Dec 20 '14 at 4:42
  • $\begingroup$ @ Rahul, Thank you for your nice example. I just acquire the use of Contours->5 :). BTW, what is your mean by "based on which side of the line has the larger function value"? As we known, for streamline or the contour lines in ContourPlot, the value of each curve is a constant. $\endgroup$ – Enter Dec 20 '14 at 4:48
5
$\begingroup$

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

Second plot with arrows heads aligned

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]]]})
$\endgroup$
  • $\begingroup$ Could you pls tell me what is the meaning of ` 1, 2, 1, 3, 1, 4, 1` in Reverse[sl[[1, 2, 1, 3, 1, 4, 1]]], especially, the last 1. Thank you! $\endgroup$ – Enter Dec 22 '14 at 3:18
  • $\begingroup$ 1, 2, 1, 3, 1, 4 is the Position of Line[x_], and the extra 1 at the end gives the position of x_ itself, which is what we want to Reverse. However, it is better to use the expression in Addendum, which is more flexible. $\endgroup$ – bbgodfrey Dec 22 '14 at 3:41
  • $\begingroup$ @ bbgodfrey. Yes, Thank you! Let me take some time to understand it. By the way, your caution about the accuracy of hxxx gives me a idea to eliminate the 3rd derivative in my function by manipulate my original equation. Thus, this problem should be solved by a mathematical trick which is out of the range of ContourPlot. $\endgroup$ – Enter Dec 22 '14 at 3:51
6
$\begingroup$

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]}

Mathematica graphics

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:

Mathematica graphics

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

$\endgroup$
  • $\begingroup$ @ Michael E2, Thank you very much!! $\endgroup$ – Enter Dec 20 '14 at 12:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.