Index: /issm/trunk/src/m/classes/public/flowlines.m
===================================================================
--- /issm/trunk/src/m/classes/public/flowlines.m	(revision 2561)
+++ /issm/trunk/src/m/classes/public/flowlines.m	(revision 2561)
@@ -0,0 +1,155 @@
+function flowpath=flowlines(index,x,y,u,v,x0,y0,varargin)
+%FLOWLINES - compute flowlines from a given set of seed points
+%
+%   Usage:
+%      flowpath=flowlines(index,x,y,u,v,x0,y0)
+%
+%   the velocity field is given by the couple (u,v) and the coordinates
+%   of the seed points are (x0,y0). One can use one or several seed 
+%   points
+%
+%   Example:
+%      flowpath=flowlines(md.elements,md.x,md.y,md.vx,md.vy,x0,y0)
+
+%check input size
+if nargin>8 | nargin<7,
+	help flowlines
+	error('flowlines error message: bad usage');
+end
+
+%check input
+if (length(x)~=length(y) | length(x)~=length(u) | length(x)~=length(v)),
+	error('flowlines error message: x,y,u and v must have the same length');
+end
+if length(x)<3,
+	error('flowlines error message: at least one element is required');
+end
+if length(x0)~=length(y0),
+	error('flowlines error message: x0 and y0 do not have the same length');
+end
+
+%get maxiter and precision
+if nargin==8
+	maxiter=varargin{1};
+else
+	maxiter=200; %maximum number of iterations
+end
+precision=1; %division of each segment (higer precision increases number of segments)
+
+%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('x',cell(N,1),'y',cell(N,1));
+for i=1:N,
+	flowpath(i).x=x0(i);
+	flowpath(i).y=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 ));
+
+%take velocity for each element
+u=u(index)*[1;1;1]/3;
+v=v(index)*[1;1;1]/3;
+
+%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))).x(end)=[];
+		flowpath(queue(listnan(i))).y(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)).x(end)+ut(i)*length_tria(tria(i))/precision;
+		Y(queue(i))=flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision;
+		flowpath(queue(i)).x=[flowpath(queue(i)).x;flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision];
+		flowpath(queue(i)).y=[flowpath(queue(i)).y;flowpath(queue(i)).y(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))).x(1)=[];
+		flowpath(queue(listnan(i))).y(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)).x(1)+ut(i)*length_tria(tria(i))/precision;
+		Y(queue(i))=flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision;
+		flowpath(queue(i)).x=[flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).x];
+		flowpath(queue(i)).y=[flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).y];
+	end
+end
Index: /issm/trunk/src/m/classes/public/plot/plot_streamlines.m
===================================================================
--- /issm/trunk/src/m/classes/public/plot/plot_streamlines.m	(revision 2560)
+++ /issm/trunk/src/m/classes/public/plot/plot_streamlines.m	(revision 2561)
@@ -4,8 +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,8 +15,4 @@
 	return
 end
-
-%take velocity for each element
-u=u(index)*[1;1;1]/3;
-v=v(index)*[1;1;1]/3;
 
 %initialize flowpath
@@ -37,121 +29,10 @@
 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