4
$\begingroup$

I have a Navier-Stokes solution for a flow past a cylinder. I have attached the DAT file. The columns in the DAT file contain the x,y,u,v,p data, respectively. The DAT file can be obtained from our attachment repository here.

I am trying to make Mathematica plot a contour plot that looks like this for velocity and pressure

Velocity field

Pressure

But when I use the following to plot my u,v,p fields, the contours erase my geometry. This is what I have done so far:

Data = Import[
   "https://raw.github.com/stackmma/Attachments/master/data_2571_187.txt",
   "Table", "IgnoreEmptyLines" -> True];
x = Data[[All, 1]];
y = Data[[All, 2]];
u = Data[[All, 3]];
v = Data[[All, 4]];
p =  Data[[All, 5]];
gridsolution = Transpose[{x, y}];
velocityu = Transpose[{x, y, u}];
velocityv = Transpose[{x, y, v}];
pressure = Transpose[{x, y, p}];
ListContourPlot[velocityu, PlotRange -> All, Contours -> 99, 
ContourStyle -> {None}, ColorFunction -> "Rainbow", 
PlotLegends -> Automatic, PlotLabel -> "u(x,y)", 
AxesLabel -> {"x", "y"}]
ListContourPlot[velocityv, PlotRange -> All, Contours -> 99, 
ContourStyle -> {None}, ColorFunction -> "Rainbow", 
PlotLegends -> Automatic, PlotLabel -> "v(x,y)", 
AxesLabel -> {"x", "y"}]
ListContourPlot[pressure, PlotRange -> All, Contours -> 99, 
ContourStyle -> {None}, ColorFunction -> "Rainbow", 
PlotLegends -> Automatic, PlotLabel -> "P(x,y)", 
AxesLabel -> {"x", "y"} ]

This results in the following figures:

u v pressure

How can I make my velocity and pressure plot look like the COMSOL plots above? Mathematica doesn't understand that there is a cylinder in the middle, but just overwrites the geometry with contours.

$\endgroup$
  • 4
    $\begingroup$ I can't see the .dat file. Am I missing something? $\endgroup$ – Jonathan Shock May 23 '13 at 6:32
  • $\begingroup$ if its the visual effect you are after it would be straightforward to add a white disc at the center? $\endgroup$ – chris May 23 '13 at 7:29
6
$\begingroup$

Revisited answer after your data was made available

After you made your data available an easier solution came to my mind.

pressure = Import[
  "https://raw.github.com/stackmma/Attachments/master/data_2571_187.txt",
  "Table"][[All, {1, 2, 5}]];

Graphics[{Hue[#5], Tooltip[Point[{#1, #2}], #5]} & @@@ data]

enter image description here

I included a tooltip so you can inspect the values of your data. As you see the values are positive in the regions you are interested in. Therefore, one very simple approach is to define the PlotRange accordingly:

ListContourPlot[pressure, Contours -> 20, 
 PlotRange -> {Automatic, Automatic, {10^-6, Automatic}}, Mesh -> All,
  ColorFunction -> "Rainbow"]

enter image description here

The same can be achieve even simpler with RegionFunction (I switch to ListDensityPlot as explained in the Old answer below)

ListDensityPlot[pressure, PlotRange -> All, 
 RegionFunction -> Function[{x, y, z}, z > 0], 
 ColorFunction -> "Rainbow"]

If you want a more perfect circle in the middle, you should extract the radius of the inner circle as described in the very end of this article

With[{minradius = Min[Norm /@ pressure[[All, {1, 2}]]]},
 ListDensityPlot[pressure, PlotRange -> All, 
  RegionFunction -> Function[{x, y, z}, Norm[{x, y}] > minradius], 
  ColorFunction -> "Rainbow"]
 ]

Therefore, opposed to what I explain in the next section, for your specific data there is no need for an Interpolation. I'm sorry that this answer is such a mess but without your data in the first place it was more guessing.

Old answer

To explain what happens let's consider a simple example. I'm creating a list of tuples $\{\{x_1,y_1,p_1\},\{x_2,y_2,p_2\},\ldots\}$ randomly distributed but sparing out the inner circle like your data:

data = Last@Last@Reap@Table[
      With[{p = RandomReal[{-2, 2}, 2]},
       If[Norm[p] < 1 && Norm[p] > .2, 
        Sow@Append[p, Sin[Times @@ p]]]], {10000}];

Graphics[{Hue[#3], Point[{#1, #2}]} & @@@ data]

enter image description here

I don't know whether your data is on a regular grid but that's not the important thing right now.

First of all, what you want is probably not a contour plot but a density plot. Let us try a simple ListDensityPlot using Mesh->All to see what happens:

ListDensityPlot[data, ColorFunction -> Hue, 
 ColorFunctionScaling -> False, Mesh -> All]

enter image description here

As you see, in order to create a colorized area Mathematica has triangulated your data points. Now here comes the clue: How should Mathematica know that you want to leave out the center? It just thinks that you probably have not many measured points there and therefore it does, what it is doing for the rest of the area: It interpolates. And since the center is surrounded by measured points, this is pretty easy. (This is opposed to outside the circle, where interpolation is not possible!)

To solve your problem, I would suggest you use Interpolation to interpolate your data, use DensityPlot to draw it and the option RegionFunction to exclude whatever region you don't want in your plot:

With[{ip = Interpolation[data, InterpolationOrder -> 1]},
 func[p : {x_, y_}] := If[Norm[p] < 1 && Norm[p] > .2, ip[x, y], 0]
 ]

DensityPlot[func[{x, y}], {x, -1, 1}, {y, -1, 1}, 
 ColorFunction -> Hue, ColorFunctionScaling -> False, 
 RegionFunction -> 
  Function[{x, y, z}, Norm[{x, y}] < .98 && Norm[{x, y}] > .2],
 PlotPoints -> 50, MaxRecursion -> 1]

enter image description here

Note that ListDensityPlot has the RegionFunction option too, but there is a pitfall hidden. ListContourPlot seems to use not really a region. Instead it just does not take data points outside outside the region into account. This gives surprising results. Try yourself.

Side note: If you don't know the values of your inner circle, you can just calculate it by using Norm and Min on all your {x,y} pairs.

$\endgroup$

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.