利用matlab绘制二维均匀流线和向量场
利用matlab繪制二維均勻流線和向量場(chǎng)(向量場(chǎng)彩色箭頭,顏色隨變量變化)
- 0前言
- 1 均勻流線的繪制
- 2 繪制彩色的短線圖
- 3 繪制彩色的均勻流線
- 4 運(yùn)動(dòng)的彩色箭頭流線圖
0前言
之前一篇文章matlab流場(chǎng)可視化后處理我簡(jiǎn)單介紹了很多繪制流場(chǎng)可視化的方法,但是在流線方面一直不是很滿意,因?yàn)閙atlab自帶的繪制流線函數(shù)streamslice沒有開源,有些參數(shù)(比如每根流線的顏色)不能修改等,所以我嘗試著自己寫一種繪制二維均勻流線的程序,以便滿足更多的需求。
其中,matlab官方file exchange里,也有一些大神的代碼給出了很好的結(jié)果:
Evenly Spaced Streamlines
Streamcolor
StreakArrow
關(guān)于如何用matlab繪制均勻的流線,我采用的方法是下面這篇論文提到的算法:
Jobard B., Lefer W. (1997) Creating Evenly-Spaced Streamlines of Arbitrary Density. In: Lefer W., Grave M. (eds) Visualization in Scientific Computing ’97. Eurographics. Springer, Vienna
1 均勻流線的繪制
論文Creating Evenly-Spaced Streamlines of Arbitrary Density中提出的繪制流線的算法可以描述為:
先初始化一根流線,然后在這個(gè)流線周圍逐漸的擴(kuò)展新的流線。文章中定義了兩個(gè)距離,一個(gè)是用于產(chǎn)生新流線種子點(diǎn)的距離dstart,一個(gè)是用于結(jié)束流線生成的距離dend。新的流線擴(kuò)展方式為,在距離舊的流線為dstart的隨機(jī)位置上,生成一個(gè)種子點(diǎn),對(duì)該種子點(diǎn)進(jìn)行新流線的生成(同時(shí)向頭和尾段生成流線),當(dāng)新流線兩端的點(diǎn)距某個(gè)舊流線的距離小于dend時(shí),停止生成流線,將該流線放入舊流線內(nèi)保存。當(dāng)任意位置都不能滿足新流線的生成時(shí),則認(rèn)為所有流線都已經(jīng)生成,結(jié)束計(jì)算。
用偽代碼的形式可以表述為
1 隨機(jī)繪制第一根流線,作為初始流線。把初始流線加入已生成的流線列表內(nèi)。把初始流線記為當(dāng)前流線。 2 令Finish=False,當(dāng)Finish==False時(shí),進(jìn)行循環(huán)3 選擇一個(gè)種子點(diǎn),該點(diǎn)距 當(dāng)前流線 的距離為dstart,且距離其它 已生成的流線列表內(nèi) 的流線距離大于dstart4 如果能夠找到一個(gè)符合條件的種子點(diǎn),則生成新的流線。5 新流線以種子點(diǎn)為開始,向前向后伸展。當(dāng)流線的前端新計(jì)算的流線點(diǎn)距離某條流線的距離小于dend,則停止前端流線的生成。后端同樣。6 將新流線保存到已生成的流線列表內(nèi)。7 如果不能找到符合條件的種子點(diǎn),則替換 當(dāng)前流線 ,重復(fù)48 如果 已生成的流線列表內(nèi) 所有流線都經(jīng)過步驟7嘗試一遍,則說明計(jì)算結(jié)束。令Finish==True, 停止2的循環(huán)。由于在計(jì)算中,用到了大量的距離計(jì)算,比如生成新種子點(diǎn)時(shí),需要不停的和其它所有流線計(jì)算距離,保證大于dstart,又比如在計(jì)算新流線時(shí),新流線每次生成一個(gè)點(diǎn),該點(diǎn)都需要和所有以有的流線進(jìn)行距離計(jì)算,保證距離大于dend。
這種頻繁的距離計(jì)算數(shù)量,在量級(jí)上為N^2,這導(dǎo)致當(dāng)流線上點(diǎn)比較多或者流線比較多時(shí),計(jì)算會(huì)非常緩慢。我們可以利用一些其它的方法減少這種運(yùn)算,比如利用劃分網(wǎng)格的方式,使得每次計(jì)算距離,只需要計(jì)算周圍網(wǎng)格(比如2維的周圍網(wǎng)格有8個(gè))內(nèi)所有點(diǎn)的距離,而無需計(jì)算全局所有點(diǎn)的距離?;蛘呃脴湫嗡阉魉惴?#xff0c;可以將計(jì)算量降到NlgN的量級(jí)上。
當(dāng)然,這里我利用一種更簡(jiǎn)單的方法進(jìn)行計(jì)算。利用劃分網(wǎng)格的思想,將新種子點(diǎn)的生成位置限定在網(wǎng)格上。而且,假設(shè)如果兩個(gè)點(diǎn)在同一個(gè)網(wǎng)格內(nèi),說明兩個(gè)點(diǎn)距離小于d,如果不在同一網(wǎng)格內(nèi),兩個(gè)點(diǎn)距離大于d。這種假設(shè)雖然對(duì)計(jì)算結(jié)果不精準(zhǔn),但是仍然是可接受的,而且不需要計(jì)算距離。
用偽代碼的形式可以表述為
1 生成距離dstart的網(wǎng)格。生成距離dend的網(wǎng)格。網(wǎng)格內(nèi)初始化為0 2 當(dāng)start網(wǎng)格內(nèi)所有區(qū)域都為1,則停止循環(huán)3 尋找start網(wǎng)格區(qū)域?yàn)?的區(qū)域,隨機(jī)選擇一個(gè)區(qū)域作為生成種子點(diǎn)。4 以種子點(diǎn)為起始點(diǎn),向前向后繪制流線5 查詢新生成的流線點(diǎn)對(duì)應(yīng)的end網(wǎng)格,記為當(dāng)前網(wǎng)格。6 如果當(dāng)前網(wǎng)格和上一個(gè)點(diǎn)的當(dāng)前網(wǎng)格是同一個(gè),則認(rèn)為當(dāng)前網(wǎng)格沒有移動(dòng),繼續(xù)生成下一個(gè)流線點(diǎn)7 如果當(dāng)前網(wǎng)格移動(dòng),且當(dāng)前end網(wǎng)格為1,則停止生成流線。8 如果當(dāng)前網(wǎng)格移動(dòng),且當(dāng)前end網(wǎng)格為0,則將當(dāng)前網(wǎng)格標(biāo)記為1,繼續(xù)生成下一個(gè)流線點(diǎn)9 將流線保存。流線經(jīng)過的所有start網(wǎng)格標(biāo)記為1。流線經(jīng)過的所有end網(wǎng)格標(biāo)記為1。之后重復(fù)3。matlab代碼為:
clear clc close all%載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %由于這里R2018b之前的版本會(huì)出現(xiàn)兼容性問題,所以更改 %xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%流線求解 streamline_sum=my_streamline(xx,yy,uu,vv,0.05);%繪制流線 figure(1) hold on xlim([xmin,xmax]); ylim([ymin,ymax]); xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis;for j=1:size(streamline_sum,2)if size(streamline_sum{j},1)<=1continue%忽略異常流線endplot(streamline_sum{j}(:,1),streamline_sum{j}(:,2),'k');%繪制箭頭arrow_direction=streamline_sum{j}(end,:)-streamline_sum{j}(end-1,:);plot_arrow(xy_lim,xy_ratio,streamline_sum{j}(end,:),[0,0,0],0.8,arrow_direction) end hold offfunction streamline_sum=my_streamline(x,y,u,v,dstart) %0處理前設(shè)置 %設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長(zhǎng)度) %dstart=0.05;默認(rèn)0.05 dend=0.5*dstart;%xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化,將流場(chǎng)縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin);num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;%初始化所有網(wǎng)格點(diǎn),0代表可以放置新點(diǎn),1代表已經(jīng)存在原有的點(diǎn) xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end);%1當(dāng)xy_start內(nèi)還有可放置的新點(diǎn)的位置時(shí),進(jìn)行循環(huán) k=0;%循環(huán)次數(shù),也是流線個(gè)數(shù) while ~all(xy_start,'all')k=k+1;%2隨機(jī)一個(gè)start內(nèi)網(wǎng)格點(diǎn)作為種子點(diǎn)[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));%3繪制流線streamline_i_1 = stream2(xn,yn, un, vn,x_pos_i,y_pos_i,0.2);streamline_i_2 = stream2(xn,yn,-un,-vn,x_pos_i,y_pos_i,0.2);%4以xy_end為標(biāo)準(zhǔn),刪除自相交或間隔太近的點(diǎn)。并順便標(biāo)記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1{1},xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2{1},xy_end,dend,xy_start,dstart);%5保存streamline_k=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線streamline_sum{k}=[xmin+streamline_k(:,1)*(xmax-xmin),ymin+streamline_k(:,2)*(ymax-ymin)];%從歸一化還原 end endfunction xpoint=id2axis(distance,id) %取網(wǎng)格的中點(diǎn) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 if id==1xpoint=min_distance/2; elseif id==Nxpoint=1-min_distance/2; elsexpoint=min_distance+(id-1.5)*distance; end endfunction [sl_i,xy_end,xy_start]=delete_self(sl_i,xy_end,dend,xy_start,dstart) %sl_i streamline流線,兩列N行形式 N=size(sl_i,1); pos_id_last=axis2id(sl_i(1,1),sl_i(1,2),dend); xy_end(pos_id_last)=1;%第一個(gè)點(diǎn)標(biāo)記%順便標(biāo)記xy_start pos_id_s=axis2id(sl_i(1,1),sl_i(1,2),dstart); xy_start(pos_id_s)=1; for j=2:Npos_id_now=axis2id(sl_i(j,1),sl_i(j,2),dend);if pos_id_now~=pos_id_last%如果現(xiàn)在的點(diǎn)和原有的點(diǎn)在同一區(qū)域,則不管它%如果不在同一區(qū)域,檢測(cè)新的點(diǎn)是否已經(jīng)被占用if xy_end(pos_id_now)==1%如果該點(diǎn)被占用,說明出現(xiàn)與其它流線太近的情況,則直接停止j=j-1;breakelse%如果沒被占用,則把新點(diǎn)添加上xy_end(pos_id_now)=1;pos_id_last=pos_id_now;endend%順便標(biāo)記xy_startpos_id_s=axis2id(sl_i(j,1),sl_i(j,2),dstart);xy_start(pos_id_s)=1; end sl_i(j:end,:)=[]; endfunction pos_id=axis2id(x,y,distance) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 %x的位置 if x<=min_distancepos_id_x=1; elseif x>=1-min_distancepos_id_x=N; elsepos_id_x=ceil((x-min_distance)/distance)+1; end %y的位置 if y<=min_distancepos_id_y=1; elseif y>=1-min_distancepos_id_y=N; elsepos_id_y=ceil((y-min_distance)/distance)+1; end %xy轉(zhuǎn)ind pos_id=sub2ind([N,N],pos_id_y,pos_id_x); endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀(歸一化的形狀) arrow_0=[0,0;-1,0.5;-1,-0.5]; %對(duì)方向進(jìn)行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對(duì)窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對(duì)比例尺歸一化 arrow_3=arrow_2.*xy_ratio_n+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') end生成的流線圖為:
2 繪制彩色的短線圖
這個(gè)繪制思路很簡(jiǎn)單,就是取每一個(gè)點(diǎn),然后繪制短流線即可。
代碼如下:
繪制出的效果如下:
這些流線基本都是長(zhǎng)短相同的,因?yàn)閟tream2函數(shù)生成的流線,流線上的點(diǎn)基本是等間距的。如果想做出速度V大的流線箭頭更長(zhǎng),速度小的流線更短的效果,這里用自己編寫的二階RK流線生成算法,進(jìn)行流線的繪制。由于interp2速度太慢,這里用griddedInterpolant函數(shù)進(jìn)行替換,優(yōu)化速度將近10倍(我也沒搞清楚為什么interp2速度這么慢)。
matlab代碼如下:
clear clc close all%繪制彩色短線圖(長(zhǎng)短不同) %載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %xmin=min(xx,[],'all');xmax=max(xx,[],'all'); %ymin=min(yy,[],'all');ymax=max(yy,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));V2=sqrt(uu.^2+vv.^2);%計(jì)算速度 %繪制變量顏色條 N_color=32; P=V2;%指定變量為V2 %P_max=max(P,[],'all');P_min=min(P,[],'all'); P_max=max(max(max(P)));P_min=min(min(min(P)));Num = numel(xx);num_streamline=20;%每條流線上的點(diǎn)數(shù)量 V2_max=max(max(max(V2))); dt=3.0*min([xx(1,2)-xx(1,1),yy(2,1)-yy(1,1)])/V2_max/num_streamline;for k=1:Numstartx=xx(k);starty=yy(k);sl_i=stream2_RK2(xx,yy,uu,vv,startx,starty,dt,num_streamline); sl_sum{k}=sl_i; endmcp=colormap(parula(N_color)); [~,~,P_color] = histcounts(P(:),linspace(P_min,P_max,N_color)); mlw=linspace(0.5,2,N_color);figure(1) hold on xlim([xmin,xmax]) ylim([ymin,ymax]) xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis;for k=1:Numplot(sl_sum{k}(:,1),sl_sum{k}(:,2),'color',mcp(P_color(k),:),'linewidth',mlw(P_color(k)))if size(sl_sum{k},1)<=1continueend%繪制箭頭arrow_direction=sl_sum{k}(end,:)-sl_sum{k}(end-1,:);plot_arrow(xy_lim,xy_ratio,sl_sum{k}(end,:),mcp(P_color(k),:),0.5*mlw(P_color(k)),arrow_direction) end hold off caxis([P_min,P_max]) colorbarfunction streamline_i=stream2_RK2(x,y,u,v,startx,starty,dt,N) streamline_i=zeros(N,2); streamline_i(1,:)=[startx,starty]; x_old=startx;y_old=starty; F_u = griddedInterpolant(x',y',u','linear'); F_v = griddedInterpolant(x',y',v','linear'); for k=2:N%利用改進(jìn)歐拉法(或者叫2階Runge-Kutta,預(yù)估校正)%interp2太慢,放棄 % u_K1=interp2(x,y,u,x_old,y_old,'linear')*dt; % v_K1=interp2(x,y,v,x_old,y_old,'linear')*dt;u_K1 = F_u(x_old,y_old)*dt;v_K1 = F_v(x_old,y_old)*dt; % u_K2=interp2(x,y,u,x_old+0.5*u_K1,y_old+0.5*v_K1,'linear')*dt; % v_K2=interp2(x,y,v,x_old+0.5*u_K1,y_old+0.5*v_K1,'linear')*dt;u_K2 = F_u(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;v_K2 = F_v(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;x_new=x_old+0.5*(u_K1+u_K2);y_new=y_old+0.5*(v_K1+v_K2);%保存streamline_i(k,:)=[x_new,y_new];x_old=x_new;y_old=y_new;if isnan(x_new) || isnan(y_new)streamline_i(k+1:end,:)=[];breakend end endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀(歸一化的形狀) arrow_0=[0,0;-1,0.5;-1,-0.5]; %對(duì)方向進(jìn)行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對(duì)窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對(duì)比例尺歸一化 arrow_3=arrow_2.*xy_ratio_n+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') end生成的流線效果如下:
3 繪制彩色的均勻流線
流線的繪制采用第1小節(jié)方法,但是采用控制stream2函數(shù)的點(diǎn)數(shù)來控制流線的長(zhǎng)度(沒有采用第2小節(jié)的RK方法)。顏色的繪制采用第2小節(jié)的方法。
代碼如下:
clear clc close all%繪制彩色長(zhǎng)線圖(長(zhǎng)短不同) %載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %xmin=min(xx,[],'all');xmax=max(xx,[],'all'); %ymin=min(yy,[],'all');ymax=max(yy,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%流線求解 [streamline_sum,streamline_seed]=my_streamline_mutli(xx,yy,uu,vv,0.03,32);%繪制變量顏色條 V2=sqrt(uu.^2+vv.^2);N_color=16;%流線的顏色種類 P=V2;%把速度作為顏色的變量 %P_max=max(P,[],'all');P_min=min(P,[],'all'); P_max=max(max(max(P)));P_min=min(min(min(P)));mcp=colormap(jet(N_color)); P_seed=interp2(xx,yy,P,streamline_seed(:,1),streamline_seed(:,2));[~,~,P_color] = histcounts(P_seed(:),linspace(P_min,P_max,N_color)); mlw=linspace(0.5,2,N_color);%繪制流線 figure(1) hold on for j=1:size(streamline_sum,2)plot(streamline_sum{j}(:,1),streamline_sum{j}(:,2),'color',mcp(P_color(j),:)...,'linewidth',mlw(P_color(j)));%繪制流線 end xlim([xmin,xmax]) ylim([ymin,ymax]) caxis([P_min,P_max]) colorbar %繪制箭頭 xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis; for k=1:size(streamline_sum,2)%繪制箭頭arrow_direction=streamline_sum{k}(end,:)-streamline_sum{k}(end-1,:);plot_arrow(xy_lim,xy_ratio,streamline_sum{k}(end,:),mcp(P_color(k),:),mlw(P_color(k)),arrow_direction) end hold offfunction [streamline_sum,streamline_seed]=my_streamline_mutli(x,y,u,v,dstart,num) %0處理前設(shè)置 %設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長(zhǎng)度) %dstart=0.05;默認(rèn)0.05 dend=0.5*dstart;%xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化,將流場(chǎng)縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin);num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;%初始化所有網(wǎng)格點(diǎn),0代表可以放置新點(diǎn),1代表已經(jīng)存在原有的點(diǎn) xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end);%將流線劃分為num種,速度越大的流線越長(zhǎng) length_sl=linspace(5,40,num); V2=sqrt(un.^2+vn.^2); V2_max=max(max(max(V2))); V2_min=min(min(min(V2))); %V2_min=min(V2,[],'all');V2_max=max(V2,[],'all');V2_space=linspace(V2_min,V2_max,num+1);%1當(dāng)xy_start內(nèi)還有可放置的新點(diǎn)的位置時(shí),進(jìn)行循環(huán) k=0;%循環(huán)次數(shù),也是流線個(gè)數(shù) while ~all(xy_start,'all')k=k+1;%2隨機(jī)一個(gè)start內(nèi)網(wǎng)格點(diǎn)作為種子點(diǎn)[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));streamline_seed(k,:)=[x_pos_i,y_pos_i];%保存種子點(diǎn)V2_seed=interp2(xn,yn,V2,x_pos_i,y_pos_i);%計(jì)算種子點(diǎn)處的速度[~,~,sl_N] = histcounts(V2_seed,V2_space);num_streamline=round(length_sl(sl_N));%3繪制流線streamline_i_1 = stream2(xn,yn, un, vn,x_pos_i,y_pos_i,[0.1,num_streamline]);streamline_i_2 = stream2(xn,yn,-un,-vn,x_pos_i,y_pos_i,[0.1,num_streamline]);%4以xy_end為標(biāo)準(zhǔn),刪除自相交或間隔太近的點(diǎn)。并順便標(biāo)記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1{1},xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2{1},xy_end,dend,xy_start,dstart);%5保存streamline_k=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線streamline_sum{k}=[xmin+streamline_k(:,1)*(xmax-xmin),ymin+streamline_k(:,2)*(ymax-ymin)];%從歸一化還原 end streamline_seed=[streamline_seed(:,1)*(xmax-xmin)+xmin,streamline_seed(:,2)*(ymax-ymin)+ymin]; endfunction xpoint=id2axis(distance,id) %取網(wǎng)格的中點(diǎn)N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 if id==1xpoint=min_distance/2; elseif id==Nxpoint=1-min_distance/2; elsexpoint=min_distance+(id-1.5)*distance; end endfunction [sl_i,xy_end,xy_start]=delete_self(sl_i,xy_end,dend,xy_start,dstart) %sl_i streamline流線,兩列N行形式 N=size(sl_i,1); pos_id_last=axis2id(sl_i(1,1),sl_i(1,2),dend); xy_end(pos_id_last)=1;%第一個(gè)點(diǎn)標(biāo)記%順便標(biāo)記xy_start pos_id_s=axis2id(sl_i(1,1),sl_i(1,2),dstart); xy_start(pos_id_s)=1; for j=2:Npos_id_now=axis2id(sl_i(j,1),sl_i(j,2),dend);if pos_id_now~=pos_id_last%如果現(xiàn)在的點(diǎn)和原有的點(diǎn)在同一區(qū)域,則不管它%如果不在同一區(qū)域,檢測(cè)新的點(diǎn)是否已經(jīng)被占用if xy_end(pos_id_now)==1%如果該點(diǎn)被占用,說明出現(xiàn)與其它流線太近的情況,則直接停止j=j-1;breakelse%如果沒被占用,則把新點(diǎn)添加上xy_end(pos_id_now)=1;pos_id_last=pos_id_now;endend%順便標(biāo)記xy_startpos_id_s=axis2id(sl_i(j,1),sl_i(j,2),dstart);xy_start(pos_id_s)=1; end sl_i(j:end,:)=[]; endfunction pos_id=axis2id(x,y,distance) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離%x的位置 if x<=min_distancepos_id_x=1; elseif x>=1-min_distancepos_id_x=N; elsepos_id_x=ceil((x-min_distance)/distance)+1; end%y的位置 if y<=min_distancepos_id_y=1; elseif y>=1-min_distancepos_id_y=N; elsepos_id_y=ceil((y-min_distance)/distance)+1; end%xy轉(zhuǎn)ind pos_id=sub2ind([N,N],pos_id_y,pos_id_x); endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀(歸一化的形狀) arrow_0=[0,0;-1,0.5;-1,-0.5]; %對(duì)方向進(jìn)行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對(duì)窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對(duì)比例尺歸一化 %arrow_3=arrow_2.*xy_ratio_n+xy_arrow;%由于兼容性問題,更改為下面語句 xy_ratio_n2=ones(3,1)*xy_ratio_n; arrow_3=arrow_2.*xy_ratio_n2+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') end4 運(yùn)動(dòng)的彩色箭頭流線圖
一開始我本來是打算用第2小節(jié)的方式繪制一個(gè)動(dòng)圖,但是發(fā)現(xiàn)原來的箭頭(粒子)離開之后,沒有新的箭頭進(jìn)行補(bǔ)充。如果明確流場(chǎng)的入口和出口,可以采用出口離開多少粒子,就在入口添加多少粒子的方式,維持流場(chǎng)內(nèi)數(shù)量恒定。但是對(duì)于復(fù)雜流場(chǎng)則不容易做到這點(diǎn)。
因此,我采用了均勻流線的方式,每次運(yùn)動(dòng)完成后,刪除離得太近的流線,然后在空白處補(bǔ)充新的流線。效果如下:
這樣做的優(yōu)點(diǎn)在于每一幀都滿足均勻流線的要求,而且保證絕大多數(shù)流線不會(huì)憑空的產(chǎn)生和消失,即使是消失也是逐漸消失。
代碼如下:
clear clc close all %繪制彩色長(zhǎng)線圖(長(zhǎng)短不同,而且是會(huì)動(dòng)的那種gif圖)%載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %xmin=min(xx,[],'all');xmax=max(xx,[],'all'); %ymin=min(yy,[],'all');ymax=max(yy,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));dstart=0.03;%網(wǎng)格寬度%流線求解 [streamline_sum,streamline_seed]=my_streamline_mutli(xx,yy,uu,vv,dstart,32);%繪制變量顏色條 V2=sqrt(uu.^2+vv.^2);N_color=16;%流線的顏色種類 P=V2;%把速度作為顏色的變量 %P_max=max(P,[],'all');P_min=min(P,[],'all'); P_max=max(max(max(P)));P_min=min(min(min(P)));mcp=colormap(jet(N_color)); P_seed=interp2(xx,yy,P,streamline_seed(:,1),streamline_seed(:,2));[~,~,P_color] = histcounts(P_seed(:),linspace(P_min,P_max,N_color)); mlw=linspace(0.5,2,N_color);%繪制流線 figure(1) hold on for j=1:size(streamline_sum,2)plot(streamline_sum{j}(:,1),streamline_sum{j}(:,2),'color',mcp(P_color(j),:)...,'linewidth',mlw(P_color(j)));%繪制流線 end xlim([xmin,xmax]) ylim([ymin,ymax])%固定圖窗大小 caxis([P_min,P_max])%設(shè)定顏色條范圍 colorbar%顯示色條xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis; for k=1:size(streamline_sum,2)%繪制箭頭if size(streamline_sum{k},1)<=2continueendarrow_direction=streamline_sum{k}(end,:)-streamline_sum{k}(end-1,:);plot_arrow(xy_lim,xy_ratio,streamline_sum{k}(end,:),mcp(P_color(k),:),mlw(P_color(k)),arrow_direction) end hold off%保存gif圖 frame=getframe(gca); im=frame2im(frame); [I,map]=rgb2ind(im,36); imwrite(I,map,'temp.gif','gif', 'Loopcount',inf,'DelayTime',0.08);%第一張%后續(xù)的變化 for k=1:60%延長(zhǎng)dt秒之后的圖像streamline_sum=my_streamline_time(xx,yy,uu,vv,streamline_sum,dstart,0.05);figure(1)clfhold onfor j=1:size(streamline_sum,2)sl_i=streamline_sum{j};if isempty(sl_i)continueendL_sl_i=size(sl_i,1);P_sl=interp2(xx,yy,P,sl_i(round(L_sl_i/2),1),sl_i(round(L_sl_i/2),2));[~,~,P_color] = histcounts(P_sl,linspace(P_min,P_max,N_color));P_color(P_color==0)=1;plot(sl_i(:,1),sl_i(:,2),'color',mcp(P_color,:)...,'linewidth',mlw(P_color));%繪制流線%繪制箭頭if size(sl_i,1)<=2continueendarrow_direction=sl_i(end,:)-sl_i(end-1,:);plot_arrow(xy_lim,xy_ratio,sl_i(end,:),mcp(P_color,:),mlw(P_color),arrow_direction)endhold offxlim([xmin,xmax])ylim([ymin,ymax])%固定圖窗大小caxis([P_min,P_max])%設(shè)定顏色條范圍colorbar%顯示色條pause(0.5)%保存gif圖frame=getframe(gca);im=frame2im(frame);[I,map]=rgb2ind(im,36);imwrite(I,map,'temp.gif','gif','WriteMode','append','DelayTime',0.08); endfunction [streamline_sum,streamline_seed]=my_streamline_mutli(x,y,u,v,dstart,num) %繪制流線 %0處理前設(shè)置%設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長(zhǎng)度) %dstart=0.05;默認(rèn)0.05 dend=0.5*dstart; %xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化,將流場(chǎng)縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin); F_u = griddedInterpolant(xn',yn',un','linear'); F_v = griddedInterpolant(xn',yn',vn','linear'); F_u_n = griddedInterpolant(xn',yn',-un','linear'); F_v_n = griddedInterpolant(xn',yn',-vn','linear'); num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;%初始化所有網(wǎng)格點(diǎn),0代表可以放置新點(diǎn),1代表已經(jīng)存在原有的點(diǎn) xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end); %將流線劃分為num種,速度越大的流線越長(zhǎng) length_sl=linspace(5,40,num); V2=sqrt(un.^2+vn.^2); %V2_min=min(V2,[],'all');V2_max=max(V2,[],'all'); V2_max=max(max(max(V2))); V2_min=min(min(min(V2)));V2_space=linspace(V2_min,V2_max,num+1);%1當(dāng)xy_start內(nèi)還有可放置的新點(diǎn)的位置時(shí),進(jìn)行循環(huán) k=0;%循環(huán)次數(shù),也是流線個(gè)數(shù) while ~all(xy_start,'all')k=k+1;%2隨機(jī)一個(gè)start內(nèi)網(wǎng)格點(diǎn)作為種子點(diǎn)[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));streamline_seed(k,:)=[x_pos_i,y_pos_i];%保存種子點(diǎn)%3繪制流線streamline_i_1=stream2_RK2(F_u ,F_v ,x_pos_i,y_pos_i,0.01,20);streamline_i_2=stream2_RK2(F_u_n,F_v_n,x_pos_i,y_pos_i,0.01,20);%4以xy_end為標(biāo)準(zhǔn),刪除自相交或間隔太近的點(diǎn)。并順便標(biāo)記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1,xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2,xy_end,dend,xy_start,dstart);%5保存streamline_k=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線streamline_sum{k}=[xmin+streamline_k(:,1)*(xmax-xmin),ymin+streamline_k(:,2)*(ymax-ymin)];%從歸一化還原 end streamline_seed=[streamline_seed(:,1)*(xmax-xmin)+xmin,streamline_seed(:,2)*(ymax-ymin)+ymin]; endfunction xpoint=id2axis(distance,id) %取網(wǎng)格的中點(diǎn) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 if id==1xpoint=min_distance/2; elseif id==Nxpoint=1-min_distance/2; elsexpoint=min_distance+(id-1.5)*distance; end endfunction [sl_i,xy_end,xy_start]=delete_self(sl_i,xy_end,dend,xy_start,dstart) %sl_i streamline流線,兩列N行形式 N=size(sl_i,1); pos_id_last=axis2id(sl_i(1,1),sl_i(1,2),dend); xy_end(pos_id_last)=1;%第一個(gè)點(diǎn)標(biāo)記%順便標(biāo)記xy_start pos_id_s=axis2id(sl_i(1,1),sl_i(1,2),dstart); xy_start(pos_id_s)=1; for j=2:Npos_id_now=axis2id(sl_i(j,1),sl_i(j,2),dend);if pos_id_now~=pos_id_last%如果現(xiàn)在的點(diǎn)和原有的點(diǎn)在同一區(qū)域,則不管它%如果不在同一區(qū)域,檢測(cè)新的點(diǎn)是否已經(jīng)被占用if xy_end(pos_id_now)==1%如果該點(diǎn)被占用,說明出現(xiàn)與其它流線太近的情況,則直接停止j=j-1;breakelse%如果沒被占用,則把新點(diǎn)添加上xy_end(pos_id_now)=1;pos_id_last=pos_id_now;endend%順便標(biāo)記xy_startpos_id_s=axis2id(sl_i(j,1),sl_i(j,2),dstart);xy_start(pos_id_s)=1; end sl_i(j:end,:)=[];%刪除停止之后剩余的流線 endfunction pos_id=axis2id(x,y,distance) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 %x的位置 if x<=min_distancepos_id_x=1; elseif x>=1-min_distancepos_id_x=N; elsepos_id_x=ceil((x-min_distance)/distance)+1; end %y的位置 if y<=min_distancepos_id_y=1; elseif y>=1-min_distancepos_id_y=N; elsepos_id_y=ceil((y-min_distance)/distance)+1; end %xy轉(zhuǎn)ind pos_id=sub2ind([N,N],pos_id_y,pos_id_x); endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀(歸一化的形狀) arrow_0=[0,0;-1,0.5;-1,-0.5]; %對(duì)方向進(jìn)行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對(duì)窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對(duì)比例尺歸一化 %arrow_3=arrow_2.*xy_ratio_n+xy_arrow;%由于兼容性問題,更改為下面語句 xy_ratio_n2=ones(3,1)*xy_ratio_n; arrow_3=arrow_2.*xy_ratio_n2+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') endfunction streamline_i=stream2_RK2(F_u,F_v,startx,starty,dt,N) %繪制流線,采用RK方法 %利用改進(jìn)歐拉法(或者叫2階Runge-Kutta,預(yù)估校正) streamline_i=zeros(N,2); streamline_i(1,:)=[startx,starty]; x_old=startx;y_old=starty; for k=2:Nu_K1 = F_u(x_old,y_old)*dt;v_K1 = F_v(x_old,y_old)*dt;u_K2 = F_u(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;v_K2 = F_v(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;x_new=x_old+0.5*(u_K1+u_K2);y_new=y_old+0.5*(v_K1+v_K2);%保存streamline_i(k,:)=[x_new,y_new];x_old=x_new;y_old=y_new;if isnan(x_new) || isnan(y_new)streamline_i(k+1:end,:)=[];breakendif 0>x_new || 0>y_new || x_new>1 || y_new>1%由于插值結(jié)果都是在0-1區(qū)間內(nèi)的歸一化數(shù)據(jù),所以超出邊界的刪除。streamline_i(k+1:end,:)=[];breakend end endfunction streamline_sum=my_streamline_time(x,y,u,v,streamline_sum,dstart,dt) %0處理前設(shè)置 %xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化,將流場(chǎng)縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin);streamline_sum(cellfun(@isempty,streamline_sum))=[];%刪除空數(shù)組 %排序,優(yōu)先保證長(zhǎng)線段不被刪除 L_sl=zeros(length(streamline_sum),1); for k=1:length(streamline_sum)L_sl(k)=sum((streamline_sum{k}(end,:)-streamline_sum{k}(1,:)).^2); end [~,I]=sort(L_sl,'descend' ); streamline_sum(:,[1:length(streamline_sum)])=streamline_sum(:,[I]);%設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長(zhǎng)度) %dstart=0.05;默認(rèn)0.05 dend=0.5*dstart;%每一條流線隨時(shí)間的變化 F_u = griddedInterpolant(xn',yn',un','linear'); F_v = griddedInterpolant(xn',yn',vn','linear'); F_u_n = griddedInterpolant(xn',yn',-un','linear'); F_v_n = griddedInterpolant(xn',yn',-vn','linear');num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;for k=1:length(streamline_sum)sl_i=streamline_sum{k};if isempty(sl_i)continueendif any(isnan(sl_i(end,:)))sl_i(1,:)=[];elsesl_i(1,:)=[];if isempty(sl_i)continueendsl_i_n=(sl_i-[xmin,ymin])./[(xmax-xmin),(ymax-ymin)];%歸一化流線sl_i_now=stream2_RK2(F_u ,F_v ,sl_i_n(end,1),sl_i_n(end,2),dt,2);sl_i_n=[sl_i_n;sl_i_now(2,:)];sl_i=sl_i_n.*[(xmax-xmin),(ymax-ymin)]+[xmin,ymin];%反歸一化endstreamline_sum{k}=sl_i; end streamline_sum(cellfun(@isempty,streamline_sum))=[];%刪除空數(shù)組xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end); for k=1:length(streamline_sum)sl_i=streamline_sum{k};sl_i=flipud(sl_i);%刪除尾部,保留頭部sl_i_n=(sl_i-[xmin,ymin])./[(xmax-xmin),(ymax-ymin)];%歸一化流線%以xy_end為標(biāo)準(zhǔn),刪除自相交或間隔太近的點(diǎn)。并順便標(biāo)記xy_end[sl_i_n,xy_end,xy_start]=delete_self(sl_i_n,xy_end,dend,xy_start,dstart);sl_i=sl_i_n.*[(xmax-xmin),(ymax-ymin)]+[xmin,ymin];%反歸一化sl_i=flipud(sl_i);streamline_sum{k}=sl_i; end %之后還有一些區(qū)域太空,所以重新在那些區(qū)域生成流線 while ~all(xy_start,'all')k=k+1;%2隨機(jī)一個(gè)start內(nèi)網(wǎng)格點(diǎn)作為種子點(diǎn)[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));%3繪制流線streamline_i_1=stream2_RK2(F_u ,F_v ,x_pos_i,y_pos_i,0.01,20);streamline_i_2=stream2_RK2(F_u_n,F_v_n,x_pos_i,y_pos_i,0.01,20);%4以xy_end為標(biāo)準(zhǔn),刪除自相交或間隔太近的點(diǎn)。并順便標(biāo)記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1,xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2,xy_end,dend,xy_start,dstart);%5保存sl_i_n=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線sl_i=sl_i_n.*[(xmax-xmin),(ymax-ymin)]+[xmin,ymin];%反歸一化streamline_sum{k}=sl_i; endend總結(jié)
以上是生活随笔為你收集整理的利用matlab绘制二维均匀流线和向量场的全部?jī)?nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: jieba提取关键词时筛选词性时单词性选
- 下一篇: matlab输入二项分布函数,matla