Changeset 2561

Timestamp:

10/29/09 08:42:31 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added flowlines that compute the coordinates of a flow line that goes through a given seed point

Location:

issm/trunk/src/m/classes/public

Files:

• flowlines.m (added)
• plot/plot_streamlines.m (modified) (3 diffs)

issm/trunk/src/m/classes/public/plot/plot_streamlines.m

@@ -4 +4 @@
 %   Usage:
 %      plot_streamlines(md,options)
-
-%some options
-maxiter=200; %maximum number of iterations
-precision=1; %division of each segment (higer precision increases number of segments)
 
 %process data and model
@@ -19 +15 @@
         return
 end
-
-%take velocity for each element
-u=u(index)*[1;1;1]/3;
-v=v(index)*[1;1;1]/3;
 
 %initialize flowpath
@@ -37 +29 @@
 end
 
-%check seed points
-tria=tsearch(x,y,index,x0,y0);
-pos=find(isnan(tria));
-%seems to be a bug here: must do the test several times...
-while ~isempty(pos)
-        x0(pos)=[];
-        y0(pos)=[];
-        tria=tsearch(x,y,index,x0,y0);
-        pos=find(isnan(tria));
-end
-
-%initialize other variables
-N=length(x0);
-X=x0; Y=y0;
-flowpath=struct('xc',cell(N,1),'yc',cell(N,1),'done',cell(N,1));
-for i=1:N,
-        flowpath(i).xc=x0(i);
-        flowpath(i).yc=y0(i);
-end
-done=zeros(N,1);
-
-%get avegared length of each element
-length_tria=1/3*(sqrt( (x(index(:,1))-x(index(:,2))).^2+(y(index(:,1))-y(index(:,2))).^2 )+...
-                                                sqrt((x(index(:,1))-x(index(:,3))).^2+(y(index(:,1))-y(index(:,3))).^2 )+...
-                                                sqrt((x(index(:,2))-x(index(:,3))).^2+(y(index(:,2))-y(index(:,3))).^2 ));
-
-%initialization:
-counter=1;
-
-while any(~done) 
-
-        %find current triangle
-        queue=find(~done);
-        tria=tsearch(x,y,index,X(queue),Y(queue));
-
-        %check that the point is actually inside a triangle of the mesh
-        listnan=find(isnan(tria));
-        for i=1:length(listnan)
-                %remove the last point
-                flowpath(queue(listnan(i))).xc(end)=[];
-                flowpath(queue(listnan(i))).yc(end)=[];
-                done(queue(listnan(i)))=1;
-        end
-        tria(listnan)=[]; 
-        queue(listnan)=[];
-
-        %velocity of the current triangle and norm it
-        ut=u(tria); vt=v(tria); normv=sqrt(ut.^2+vt.^2);
-        ut=ut./normv;vt=vt./normv;
-
-        %check counter
-        if counter>maxiter
-                disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward'])
-                break
-        end
-        counter=counter+1;
-
-        %remove stagnant point
-        done(queue(find(ut==0 & vt==0)))=1;
-
-        %build next point
-        for i=1:length(queue)
-                X(queue(i))=flowpath(queue(i)).xc(end)+ut(i)*length_tria(tria(i))/precision;
-                Y(queue(i))=flowpath(queue(i)).yc(end)+vt(i)*length_tria(tria(i))/precision;
-                flowpath(queue(i)).xc=[flowpath(queue(i)).xc;flowpath(queue(i)).xc(end)+ut(i)*length_tria(tria(i))/precision];
-                flowpath(queue(i)).yc=[flowpath(queue(i)).yc;flowpath(queue(i)).yc(end)+vt(i)*length_tria(tria(i))/precision];
-        end
-end
-
-%same process but reverse (vel=-vel) to have a vcomplete flow line
-counter=1;
-X=x0; Y=y0;
-done=zeros(N,1);
-
-while any(~done) 
-
-        %find current triangle
-        queue=find(~done);
-        tria=tsearch(x,y,index,X(queue),Y(queue));
-
-        %check that the point is actually inside a triangle of the mesh
-        listnan=find(isnan(tria));
-        for i=1:length(listnan)
-                %remove the last point
-                flowpath(queue(listnan(i))).xc(1)=[];
-                flowpath(queue(listnan(i))).yc(1)=[];
-                done(queue(listnan(i)))=1;
-        end
-        tria(listnan)=[]; 
-        queue(listnan)=[];
-
-        %velocity of the current triangle and norm it
-        ut=-u(tria); vt=-v(tria); normv=sqrt(ut.^2+vt.^2);
-        ut=ut./normv;vt=vt./normv;
-
-        %check counter
-        if counter>maxiter
-                disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward'])
-                break
-        end
-        counter=counter+1;
-
-        %remove stagnant point
-        done(queue(find(ut==0 & vt==0)))=1;
-
-        %build next point
-        for i=1:length(queue)
-                X(queue(i))=flowpath(queue(i)).xc(1)+ut(i)*length_tria(tria(i))/precision;
-                Y(queue(i))=flowpath(queue(i)).yc(1)+vt(i)*length_tria(tria(i))/precision;
-                flowpath(queue(i)).xc=[flowpath(queue(i)).xc(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).xc];
-                flowpath(queue(i)).yc=[flowpath(queue(i)).yc(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).yc];
-        end
-end
+%Get flow lines
+flowpath=flowlines(index,x,y,u,v,x0,y0);
 
 %plot
 hold on
 for i=1:length(flowpath)
-        patch('Xdata',[flowpath(i).xc;NaN],'Ydata',[flowpath(i).yc;NaN],'facecolor','none','edgecolor','y');
+        patch('Xdata',[flowpath(i).x;NaN],'Ydata',[flowpath(i).y;NaN],'facecolor','none','edgecolor','y');
 end