Python equivalent of Matlab patch function

Question

bekhudi chakr 2022-3-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1666299-python-equivalent-of-matlab-patch-function

回答： Sonu Saw 2022-4-2

I am trying to code a method of characteristics code that gives a nozzle contour as well mach number patches inside the contour from Matlab to Python. I have created everything that outlines the method of characteristics and the graph plotting of the nozzle contour. The only thing I do not know how to do is plot the mach numbers' patches that is in Matlab into python. I have tried using first Scipy.interpolate grid data along with pcolormesh with nearest and linear for grid data. Nearest seemed to work well but my code will plot but the graph would take forever displaying and would just not respond so I had force quit. Linear did work but after some third iteration, I would get a QHull error, could not find what that meant. So I went to tricontourf but after a third iteration, I get the 'x and y arrays must consist of at least 3 unique points' ValueError.

Numchar is the user input for the number of characteristic lines desired. M is for mach and X and Y for x,y coordinates. Dstar is just the diameter of the throat. dataENS is for the expansion and non simple region of the nozzle while dataS is for the straightening section of the nozzle.

Matlab Code

% Starting Constant Section

% - Where the input Mach number is M1

% - Have not reached the first characteristic

ss(1,1) = 0; % \

ss(1,2) = Dstar/2; % |->

ss(1,3) = M1; % /

ss(2,1) = ss(1,1); % \

ss(2,2) = -ss(1,2); % |->

ss(2,3) = ss(1,3); % /

ss(3,1) = dataENS{1,2}.X; % \

ss(3,2) = dataENS{1,2}.Y; % |->

ss(3,3) = M1; % /

% Ending Constant Section

% - Where the output Mach number is Mexit

% - Have reached horizontal, fully expanded flow

se(1,1) = dataENS{numChar,numChar+1}.X; % \

se(1,2) = dataENS{numChar,numChar+1}.Y; % |->

se(1,3) = dataENS{numChar,numChar+1}.M; % /

se(2,1) = dataS{numChar}.X; % \

se(2,2) = 0; % |->

se(2,3) = dataS{numChar}.M; % /

se(3,1) = dataS{numChar}.X; % \

se(3,2) = dataS{numChar}.Y; % |->

se(3,3) = dataS{numChar}.M; % /

% Plot the starting and ending uniform flow patches

patch(se(:,1),se(:,2),se(:,3),'EdgeColor','none','FaceColor','interp') % Plot the upper starting patch

patch(se(:,1),-se(:,2),se(:,3),'EdgeColor','none','FaceColor','interp'); % Plot the lower starting patch

patch(ss(:,1),ss(:,2),ss(:,3),'EdgeColor','none','FaceColor','interp'); % Plot the upper ending patch

patch(ss(:,1),-ss(:,2),ss(:,3),'EdgeColor','none','FaceColor','interp'); % Plot the lower ending patch

% Expansion and non-simple region patch plotting

sol = []; % Initialize the solution array

for R = 1:1:numChar-1

for L = 1:1:numChar

if (R == L-1) % If we are on the nozzle edge

%Using tricontourf I can get as far as L = 3 in Python;

sol(1,1) = dataENS{R,L}.X;

sol(1,2) = dataENS{R,L}.Y;

sol(1,3) = dataENS{R,L}.M;

sol(2,1) = dataENS{R+1,L}.X;

sol(2,2) = dataENS{R+1,L}.Y;

sol(2,3) = dataENS{R+1,L}.M;

sol(3,1) = dataENS{R+1,L+1}.X;

sol(3,2) = dataENS{R+1,L+1}.Y;

sol(3,3) = dataENS{R+1,L+1}.M;

else % If we are not on the nozzle edge

sol(1,1) = dataENS{R,L}.X;

sol(1,2) = dataENS{R,L}.Y;

sol(1,3) = dataENS{R,L}.M;

sol(2,1) = dataENS{R+1,L}.X;

sol(2,2) = dataENS{R+1,L}.Y;

sol(2,3) = dataENS{R+1,L}.M;

sol(3,1) = dataENS{R+1,L+1}.X;

sol(3,2) = dataENS{R+1,L+1}.Y;

sol(3,3) = dataENS{R+1,L+1}.M;

sol(4,1) = dataENS{R,L+1}.X;

sol(4,2) = dataENS{R,L+1}.Y;

sol(4,3) = dataENS{R,L+1}.M;

end

patch(sol(:,1),sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

patch(sol(:,1),-sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

end

% Straightening region patch plotting

sol = [];

for L = 1:1:numChar

if (L == 1)

sol(1,1) = dataENS{numChar,L}.X;

sol(1,2) = dataENS{numChar,L}.Y;

sol(1,3) = dataENS{numChar,L}.M;

sol(2,1) = dataS{L}.X;

sol(2,2) = dataS{L}.Y;

sol(2,3) = dataS{L}.M;

sol(3,1) = dataENS{numChar,L+1}.X;

sol(3,2) = dataENS{numChar,L+1}.Y;

sol(3,3) = dataENS{numChar,L+1}.M;

else

sol(1,1) = dataENS{numChar,L}.X;

sol(1,2) = dataENS{numChar,L}.Y;

sol(1,3) = dataENS{numChar,L}.M;

sol(2,1) = dataS{L-1}.X;

sol(2,2) = dataS{L-1}.Y;

sol(2,3) = dataS{L-1}.M;

sol(3,1) = dataS{L}.X;

sol(3,2) = dataS{L}.Y;

sol(3,3) = dataS{L}.M;

sol(4,1) = dataENS{numChar,L+1}.X;

sol(4,2) = dataENS{numChar,L+1}.Y;

sol(4,3) = dataENS{numChar,L+1}.M;

end

patch(sol(:,1),sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

patch(sol(:,1),-sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

end

Here is the python equivalent I have tried

#Set up the figure

figure_title = 'Nozzle - Method of Characteristics (' + str(numChar) + ")"

fig, ax = plt.subplots()

ax.set_title(figure_title)

plt.ion

#Starting Constant Section - Where the input Mach number is M1; Have not reached the first characteristic

ss = [ [], [], [] ] #Form [ [x-coor], [y-coor], [mach nums] ]

ss[0].append(0) #\

ss[0].append(ss[0][0]) # |->

ss[0].append(dataENS_x[0,1]) #/

ss[1].append(Dstar/2) #\

ss[1].append(-ss[1][0]) # |->

ss[1].append(dataENS_y[0,1]) #/

ss[2].append(throat_mach_num) #\

ss[2].append(ss[0][2]) # |->

ss[2].append(throat_mach_num) #/

plt.tricontourf(ss[0], ss[1], ss[2],100)

ss_neg = [-ss[1][0], -ss[1][1], -ss[1][2]]

plt.tricontourf(ss[0], ss_neg, ss[2], 100)

"""

This is my attempt to use spicy.interpolate with grid data. It worked well but it would take forever to show the graph on my Mac. When I used linear, code would break in the for loop at j = 2

# Plot the starting uniform flow patches

#xis = np.linspace(min(ss[0]), max(ss[0]), 100)

#yis = np.linspace(min(ss[1]), max(ss[1]), 100)

#zis = griddata((ss[0], ss[1]), ss[2], (xis[None,:], yis[:,None]), method='nearest')

#plt.pcolormesh(xis, yis, zis, cmap='viridis')

#ss_neg_y = [-ss[1][0], -ss[1][1], -ss[1][2]] #For graphing purposes only

#yins = np.linspace(min(ss_neg_y), max(ss_neg_y), 100)

#zis2 = griddata((ss[0], ss_neg_y), ss[2], (xis[None,:], yins[:,None]), method='nearest')

#plt.pcolormesh(xis, yins, zis2, cmap='viridis') """

#Ending Constant Section - Where the output Mach number is Mexit; Have reached horizontal, fully expanded flow

se = [[], [], []] #Form [[x-coor], [y-coor], [mach nums]]

se[0].append(dataENS_x[numChar-1, numChar]) #\

se[0].append(dataS_x[numChar-1,0]) # |->

se[0].append(dataS_x[numChar-1,0]) #/

se[1].append(dataENS_y[numChar-1, numChar]) #\

se[1].append(0) # |->

se[1].append(dataS_y[numChar-1,0]) #/

se[2].append(dataENS_mach[numChar-1, numChar]) #\

se[2].append(dataS_mach[numChar-1,0]) # |->

se[2].append(dataS_mach[numChar-1,0]) #/

plt.tricontourf(se[0], se[1], se[2],100)

se_neg = [-se[1][0], -se[1][1], -se[1][2]]

plt.tricontourf(se[0], se_neg, se[2], 100)

plt.ion()

#Everything above plots well when I used tricontourf.

"""

# Plot the starting ending uniform flow patches

xise = np.linspace(min(se[0]), max(se[0]), 100)

yise = np.linspace(min(se[1]), max(se[1]), 100)

zise = griddata((se[0], se[1]), se[2], (xise[None,:], yise[:,None]), method='nearest')

plt.pcolormesh(xise, yise, zise, cmap='viridis')

se_neg_y = [-se[1][0], -se[1][1], -se[1][2]] #For graphing purposes only

yinse = np.linspace(min(se_neg_y), max(se_neg_y), 100)

zise2 = griddata((se[0], se_neg_y), se[2], (xise[None,:], yinse[:,None]), method='nearest')

plt.pcolormesh(xise, yinse, zise2, cmap='viridis')"""

#Expansion and non-simple region patch plotting

sol = [[0,0,0,0], [0,0,0,0], [0,0,0,0]] #Initialize the solution array; Form [[x-coor], [y-coor], [mach nums]]

for i in range(numChar-1):

for j in range(numChar):

if i == j-1: #If we are on the nozzle edge

sol[0][0] = dataENS_x[i,j]

sol[0][1] = dataENS_x[i+1,j]

sol[0][2] = dataENS_x[i+1,j+1]

sol[1][0] = dataENS_y[i,j]

sol[1][1] = dataENS_y[i+1,j]

sol[1][2] = dataENS_y[i+1,j+1]

sol[2][0] = dataENS_mach[i,j]

sol[2][1] = dataENS_mach[i+1,j]

sol[2][2] = dataENS_mach[i+1,j+1]

else: #If we are not on the nozzle edge

sol[0][0] = dataENS_x[i,j]

sol[0][1] = dataENS_x[i+1,j]

sol[0][2] = dataENS_x[i+1,j+1]

sol[0][3] = dataENS_x[i,j+1]

sol[1][0] = dataENS_y[i,j]

sol[1][1] = dataENS_y[i+1,j]

sol[1][2] = dataENS_y[i+1,j+1]

sol[1][3] = dataENS_y[i,j+1]

sol[2][0] = dataENS_mach[i,j]

sol[2][1] = dataENS_mach[i+1,j]

sol[2][2] = dataENS_mach[i+1,j+1]

sol[2][3] = dataENS_mach[i,j+1]

"""code would break here for both tricontourf with an error of "at least 3 unique points" and

for spicy.interpolate with grid data and pcolormesh, I would get 'Qhull' error.

Nearest seemed to work the best since the end of my code is supposed to show a ====end==== part which

dictates everything is done, which I saw, but the graph was never produced. would just load but never show and eventually freeze so I had to force quit it. Cubic was a definite no."""

plt.tricontourf(sol[0], sol[1], sol[2],100)

sol_neg = [-sol[1][0], -sol[1][1], -sol[1][2], -sol[1][3]]

plt.tricontourf(sol[0], sol_neg, sol[2],100)

plt.ioff()

"""

xisol = np.linspace(min(sol[0]), max(sol[0]), 100)

yisol = np.linspace(min(sol[1]), max(sol[1]), 100)

zisol = griddata((sol[0], sol[1]), sol[2], (xisol[None,:], yisol[:,None]), method='nearest')

plt.pcolormesh(xisol, yisol, zisol, cmap='viridis')

sol_neg_y = [-1*sol[1][0], -1* sol[1][1], -1*sol[1][2], -1*sol[1][3]]#For graphing purposes only

yinsol = np.linspace(min(sol_neg_y), max(sol_neg_y), 100)

zisol2 = griddata((sol[0], sol_neg_y), sol[2], (xisol[None,:], yinsol[:,None]), method='nearest')

plt.pcolormesh(xisol, yinsol, zisol2, cmap='viridis')"""

""" I have not reached this section since I still cannot graph all of the above section

Plotting: Patches of Mach Number

% Starting Constant Section

% - Where the input Mach number is M1

% - Have not reached the first characteristic

ss(1,1) = 0; % \

ss(1,2) = Dstar/2; % |->

ss(1,3) = M1; % /

ss(2,1) = ss(1,1); % \

ss(2,2) = -ss(1,2); % |->

ss(2,3) = ss(1,3); % /

ss(3,1) = dataENS{1,2}.X; % \

ss(3,2) = dataENS{1,2}.Y; % |->

ss(3,3) = M1; % /

% Ending Constant Section

% - Where the output Mach number is Mexit

% - Have reached horizontal, fully expanded flow

se(1,1) = dataENS{numChar,numChar+1}.X; % \

se(1,2) = dataENS{numChar,numChar+1}.Y; % |->

se(1,3) = dataENS{numChar,numChar+1}.M; % /

se(2,1) = dataS{numChar}.X; % \

se(2,2) = 0; % |->

se(2,3) = dataS{numChar}.M; % /

se(3,1) = dataS{numChar}.X; % \

se(3,2) = dataS{numChar}.Y; % |->

se(3,3) = dataS{numChar}.M; % /

% Plot the starting and ending uniform flow patches

patch(se(:,1),se(:,2),se(:,3),'EdgeColor','none','FaceColor','interp'); % Plot the upper starting patch

patch(se(:,1),-se(:,2),se(:,3),'EdgeColor','none','FaceColor','interp'); % Plot the lower starting patch

patch(ss(:,1),ss(:,2),ss(:,3),'EdgeColor','none','FaceColor','interp'); % Plot the upper ending patch

patch(ss(:,1),-ss(:,2),ss(:,3),'EdgeColor','none','FaceColor','interp'); % Plot the lower ending patch

% Expansion and non-simple region patch plotting

sol = []; % Initialize the solution array

for R = 1:1:numChar-1

for L = 1:1:numChar

if (R == L-1) % If we are on the nozzle edge

sol(1,1) = dataENS{R,L}.X;

sol(1,2) = dataENS{R,L}.Y;

sol(1,3) = dataENS{R,L}.M;

sol(2,1) = dataENS{R+1,L}.X;

sol(2,2) = dataENS{R+1,L}.Y;

sol(2,3) = dataENS{R+1,L}.M;

sol(3,1) = dataENS{R+1,L+1}.X;

sol(3,2) = dataENS{R+1,L+1}.Y;

sol(3,3) = dataENS{R+1,L+1}.M;

else % If we are not on the nozzle edge

sol(1,1) = dataENS{R,L}.X;

sol(1,2) = dataENS{R,L}.Y;

sol(1,3) = dataENS{R,L}.M;

sol(2,1) = dataENS{R+1,L}.X;

sol(2,2) = dataENS{R+1,L}.Y;

sol(2,3) = dataENS{R+1,L}.M;

sol(3,1) = dataENS{R+1,L+1}.X;

sol(3,2) = dataENS{R+1,L+1}.Y;

sol(3,3) = dataENS{R+1,L+1}.M;

sol(4,1) = dataENS{R,L+1}.X;

sol(4,2) = dataENS{R,L+1}.Y;

sol(4,3) = dataENS{R,L+1}.M;

end

patch(sol(:,1),sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

patch(sol(:,1),-sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

end

% Straightening region patch plotting

sol = [];

for L = 1:1:numChar

if (L == 1)

sol(1,1) = dataENS{numChar,L}.X;

sol(1,2) = dataENS{numChar,L}.Y;

sol(1,3) = dataENS{numChar,L}.M;

sol(2,1) = dataS{L}.X;

sol(2,2) = dataS{L}.Y;

sol(2,3) = dataS{L}.M;

sol(3,1) = dataENS{numChar,L+1}.X;

sol(3,2) = dataENS{numChar,L+1}.Y;

sol(3,3) = dataENS{numChar,L+1}.M;

else

sol(1,1) = dataENS{numChar,L}.X;

sol(1,2) = dataENS{numChar,L}.Y;

sol(1,3) = dataENS{numChar,L}.M;

sol(2,1) = dataS{L-1}.X;

sol(2,2) = dataS{L-1}.Y;

sol(2,3) = dataS{L-1}.M;

sol(3,1) = dataS{L}.X;

sol(3,2) = dataS{L}.Y;

sol(3,3) = dataS{L}.M;

sol(4,1) = dataENS{numChar,L+1}.X;

sol(4,2) = dataENS{numChar,L+1}.Y;

sol(4,3) = dataENS{numChar,L+1}.M;

end

patch(sol(:,1),sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

patch(sol(:,1),-sol(:,2),sol(:,3),'EdgeColor','none',...

'FaceColor','interp');

end

"""

The picture with the whole color scheme is the result of Matlab. That is what I want to be able to do matlab_image:https://i.stack.imgur.com/snIOz.jpg show-hairThe other picture is what I have so far with python. python_image.

My guess is that I will have to use a patches.Path where I literally do as Matlab does and draw a patch and try to find a color scheme (so if there is a way to do that, let me know).

My other guess is that I will have to import the Matlab engine into python to do the patching in which case I would need help doing that and making sure the patches are added to the python image or at least importing the image into my python graph. Or from what I read to, I can write a python package but not sure how that works or how to do it.

Besides those three guesses, I am not sure how to plot the patches. Any help is greatly appreciated.