Index: sm/trunk/src/m/classes/public/flowlines.m
===================================================================
--- /issm/trunk/src/m/classes/public/flowlines.m	(revision 5337)
+++ 	(revision )
@@ -1,160 +1,0 @@
-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),'name','','density',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
-
-%EXP compatibility
-for i=1:length(queue)
-	flowpath(queue(i)).name=['flowline' num2str(i)];
-end
Index: /issm/trunk/src/m/utils/Exp/flowlines.m
===================================================================
--- /issm/trunk/src/m/utils/Exp/flowlines.m	(revision 5338)
+++ /issm/trunk/src/m/utils/Exp/flowlines.m	(revision 5338)
@@ -0,0 +1,160 @@
+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),'name','','density',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
+
+%EXP compatibility
+for i=1:length(queue)
+	flowpath(queue(i)).name=['flowline' num2str(i)];
+end