Changeset 20535


Ignore:
Timestamp:
04/22/16 11:31:47 (9 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added more options to flowlines, that can now be used to replace downstreamflowlines.m

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  
    33%
    44%   Usage:
    5 %      flowpath=flowlines(index,x,y,u,v,x0,y0)
     5%      flowpath=flowlines(index,x,y,u,v,x0,y0,varargin)
    66%
    77%   the velocity field is given by the couple (u,v) and the coordinates
    … …  
    1111%   Example:
    1212%      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)
    1319
    1420%check input size
    … …  
    2935end
    3036
    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
     38options   = pairoptions(varargin{:});
     39maxiter   = getfieldvalue(options,'maxiter',200);
     40precision = getfieldvalue(options,'precision',1);
    3841
    3942%Create triangulation once for all and check seed points
    … …  
    6366v=v(index)*[1;1;1]/3;
    6467
    65 %initialization:
    66 counter=1;
     68if downstream,
     69        %initialization:
     70        counter=1;
    6771
    68 while any(~done)
     72        while any(~done)
    6973
    70         %find current triangle
    71         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)]);
    7377
    74         %check that the point is actually inside a triangle of the mesh
    75         listnan=find(isnan(tria));
    76         for i=1:length(listnan)
    77                 %remove the last point
    78                 flowpath(queue(listnan(i))).x(end)=[];
    79                 flowpath(queue(listnan(i))).y(end)=[];
    80                 done(queue(listnan(i)))=1;
    81         end
    82         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)=[];
    8488
    85         if isempty(tria),
    86                 break;
    87         end
     89                if isempty(tria),
     90                        break;
     91                end
    8892
    89         %velocity of the current triangle and norm it
    90         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;
    9296
    93         %check counter
    94         if counter>maxiter
    95                 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward'])
    96                 break
    97         end
    98         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;
    99103
    100         %remove stagnant point
    101         done(queue(find(ut==0 & vt==0)))=1;
     104                %remove stagnant point
     105                done(queue(find(ut==0 & vt==0)))=1;
    102106
    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
    109114        end
    110115end
    111116
    112117%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);
     118if upstream,
     119        queue=[];
     120        counter=1;
     121        X=x0; Y=y0;
     122        done=zeros(N,1);
    117123
    118 while any(~done)
     124        while any(~done)
    119125
    120         %find current triangle
    121         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)]);
    123129
    124         %check that the point is actually inside a triangle of the mesh
    125         listnan=find(isnan(tria));
    126         for i=1:length(listnan)
    127                 %remove the last point
    128                 flowpath(queue(listnan(i))).x(1)=[];
    129                 flowpath(queue(listnan(i))).y(1)=[];
    130                 done(queue(listnan(i)))=1;
    131         end
    132         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)=[];
    134140
    135         if isempty(tria),
    136                 break;
    137         end
     141                if isempty(tria),
     142                        break;
     143                end
    138144
    139         %velocity of the current triangle and norm it
    140         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;
    142148
    143         %check counter
    144         if counter>maxiter
    145                 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward'])
    146                 break
    147         end
    148         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;
    149155
    150         %remove stagnant point
    151         done(queue(find(ut==0 & vt==0)))=1;
     156                %remove stagnant point
     157                done(queue(find(ut==0 & vt==0)))=1;
    152158
    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
    159166        end
    160167end
    161168
    162 %EXP compatibility
     169%EXP compatibility (add name)
    163170for i=1:length(queue)
    164171        flowpath(queue(i)).name=['flowline' num2str(i)];
Note: See TracChangeset for help on using the changeset viewer.