Changeset 20535
- Timestamp:
- 04/22/16 11:31:47 (9 years ago)
- Location:
- issm/trunk-jpl/src/m/exp
- Files:
-
- 1 deleted
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/exp/flowlines.m
r20512 r20535 3 3 % 4 4 % Usage: 5 % flowpath=flowlines(index,x,y,u,v,x0,y0 )5 % flowpath=flowlines(index,x,y,u,v,x0,y0,varargin) 6 6 % 7 7 % the velocity field is given by the couple (u,v) and the coordinates … … 11 11 % Example: 12 12 % flowpath=flowlines(md.mesh.elements,md.mesh.x,md.mesh.y,md.vx,md.initialization.vy,x0,y0) 13 % 14 % Options: 15 % - 'maxiter': how many steps upstream and downstream of the seed points (default: 200) 16 % - 'precision': division of each segment (higer precision increases number of segments, default: 1) 17 % - 'downstream':flow line upstream of the seed points (default: 1) 18 % - 'upstream': flow line upstream of the seed points (default: 1) 13 19 14 20 %check input size … … 29 35 end 30 36 31 %get maxiter and precision 32 if nargin==8 33 maxiter=varargin{1}; 34 else 35 maxiter=200; %maximum number of iterations 36 end 37 precision=1; %division of each segment (higer precision increases number of segments) 37 %process options 38 options = pairoptions(varargin{:}); 39 maxiter = getfieldvalue(options,'maxiter',200); 40 precision = getfieldvalue(options,'precision',1); 38 41 39 42 %Create triangulation once for all and check seed points … … 63 66 v=v(index)*[1;1;1]/3; 64 67 65 %initialization: 66 counter=1; 68 if downstream, 69 %initialization: 70 counter=1; 67 71 68 while any(~done)72 while any(~done) 69 73 70 %find current triangle71 queue=find(~done);72 tria = pointLocation(trep,[X(queue),Y(queue)]);74 %find current triangle 75 queue=find(~done); 76 tria = pointLocation(trep,[X(queue),Y(queue)]); 73 77 74 %check that the point is actually inside a triangle of the mesh75 listnan=find(isnan(tria));76 for i=1:length(listnan)77 %remove the last point78 flowpath(queue(listnan(i))).x(end)=[];79 flowpath(queue(listnan(i))).y(end)=[];80 done(queue(listnan(i)))=1;81 end82 tria(listnan)=[];83 queue(listnan)=[];78 %check that the point is actually inside a triangle of the mesh 79 listnan=find(isnan(tria)); 80 for i=1:length(listnan) 81 %remove the last point 82 flowpath(queue(listnan(i))).x(end)=[]; 83 flowpath(queue(listnan(i))).y(end)=[]; 84 done(queue(listnan(i)))=1; 85 end 86 tria(listnan)=[]; 87 queue(listnan)=[]; 84 88 85 if isempty(tria),86 break;87 end89 if isempty(tria), 90 break; 91 end 88 92 89 %velocity of the current triangle and norm it90 ut=u(tria); vt=v(tria); normv=max(eps,sqrt(ut.^2+vt.^2));91 ut=ut./normv;vt=vt./normv;93 %velocity of the current triangle and norm it 94 ut=u(tria); vt=v(tria); normv=max(eps,sqrt(ut.^2+vt.^2)); 95 ut=ut./normv;vt=vt./normv; 92 96 93 %check counter94 if counter>maxiter95 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward'])96 break97 end98 counter=counter+1;97 %check counter 98 if counter>maxiter 99 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward']) 100 break 101 end 102 counter=counter+1; 99 103 100 %remove stagnant point101 done(queue(find(ut==0 & vt==0)))=1;104 %remove stagnant point 105 done(queue(find(ut==0 & vt==0)))=1; 102 106 103 %build next point 104 for i=1:length(queue) 105 X(queue(i))=flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision; 106 Y(queue(i))=flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision; 107 flowpath(queue(i)).x=[flowpath(queue(i)).x;flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision]; 108 flowpath(queue(i)).y=[flowpath(queue(i)).y;flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision]; 107 %build next point 108 for i=1:length(queue) 109 X(queue(i))=flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision; 110 Y(queue(i))=flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision; 111 flowpath(queue(i)).x=[flowpath(queue(i)).x;flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision]; 112 flowpath(queue(i)).y=[flowpath(queue(i)).y;flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision]; 113 end 109 114 end 110 115 end 111 116 112 117 %same process but reverse (vel=-vel) to have a vcomplete flow line 113 queue=[]; 114 counter=1; 115 X=x0; Y=y0; 116 done=zeros(N,1); 118 if upstream, 119 queue=[]; 120 counter=1; 121 X=x0; Y=y0; 122 done=zeros(N,1); 117 123 118 while any(~done)124 while any(~done) 119 125 120 %find current triangle121 queue=find(~done);122 tria = pointLocation(trep,[X(queue),Y(queue)]);126 %find current triangle 127 queue=find(~done); 128 tria = pointLocation(trep,[X(queue),Y(queue)]); 123 129 124 %check that the point is actually inside a triangle of the mesh125 listnan=find(isnan(tria));126 for i=1:length(listnan)127 %remove the last point128 flowpath(queue(listnan(i))).x(1)=[];129 flowpath(queue(listnan(i))).y(1)=[];130 done(queue(listnan(i)))=1;131 end132 tria(listnan)=[];133 queue(listnan)=[];130 %check that the point is actually inside a triangle of the mesh 131 listnan=find(isnan(tria)); 132 for i=1:length(listnan) 133 %remove the last point 134 flowpath(queue(listnan(i))).x(1)=[]; 135 flowpath(queue(listnan(i))).y(1)=[]; 136 done(queue(listnan(i)))=1; 137 end 138 tria(listnan)=[]; 139 queue(listnan)=[]; 134 140 135 if isempty(tria),136 break;137 end141 if isempty(tria), 142 break; 143 end 138 144 139 %velocity of the current triangle and norm it140 ut=-u(tria); vt=-v(tria); normv=max(eps,sqrt(ut.^2+vt.^2));141 ut=ut./normv;vt=vt./normv;145 %velocity of the current triangle and norm it 146 ut=-u(tria); vt=-v(tria); normv=max(eps,sqrt(ut.^2+vt.^2)); 147 ut=ut./normv;vt=vt./normv; 142 148 143 %check counter144 if counter>maxiter145 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward'])146 break147 end148 counter=counter+1;149 %check counter 150 if counter>maxiter 151 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward']) 152 break 153 end 154 counter=counter+1; 149 155 150 %remove stagnant point151 done(queue(find(ut==0 & vt==0)))=1;156 %remove stagnant point 157 done(queue(find(ut==0 & vt==0)))=1; 152 158 153 %build next point 154 for i=1:length(queue) 155 X(queue(i))=flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; 156 Y(queue(i))=flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; 157 flowpath(queue(i)).x=[flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).x]; 158 flowpath(queue(i)).y=[flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).y]; 159 %build next point 160 for i=1:length(queue) 161 X(queue(i))=flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; 162 Y(queue(i))=flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; 163 flowpath(queue(i)).x=[flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).x]; 164 flowpath(queue(i)).y=[flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).y]; 165 end 159 166 end 160 167 end 161 168 162 %EXP compatibility 169 %EXP compatibility (add name) 163 170 for i=1:length(queue) 164 171 flowpath(queue(i)).name=['flowline' num2str(i)];
Note:
See TracChangeset
for help on using the changeset viewer.