Ignore:
Timestamp:
10/30/09 10:57:03 (15 years ago)
Author:
Mathieu Morlighem
Message:

better interpolation routine call

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/classes/public/plugvelocities.m

    r2566 r2574  
    2323end
    2424
    25 %Get variables
    26 A=whos('-file',filename);
     25%load velocities
     26Names=VelFindVarNames(filename);
     27Vel=load(filename);
    2728
    28 %find x,y,vx and vy
    29 xenum=NaN; yenum=NaN; vxenum=NaN; vyenum=NaN; indexenum=NaN;
    30 if length(A)==4,
    31         isgrid=1;
    32         for i=1:4
    33                 if strcmpi(A(i).name(1),'x');
    34                         xenum=i;
    35                 elseif strcmpi(A(i).name(1),'y');
    36                         yenum=i;
    37                 else
    38                         if (strcmpi(A(i).name(end),'x') | strncmpi(A(i).name,'vx',2));
    39                                 vxenum=i;
    40                         elseif (strcmpi(A(i).name(end),'y') | strncmpi(A(i).name,'vy',2));
    41                                 vyenum=i;
    42                         end
    43                 end
    44         end
    45 elseif length(A)==5,
    46         isgrid=0;
    47         for i=1:5
    48                 if strcmpi(A(i).name(1),'x');
    49                         xenum=i;
    50                 elseif strcmpi(A(i).name(1),'y');
    51                         yenum=i;
    52                 elseif (strcmpi(A(i).name(1),'index') | strcmpi(A(i).name(1),'elements'));
    53                         indexenum=i;
    54                 else
    55                         if (strcmpi(A(i).name(end),'x') | strncmpi(A(i).name,'vx',2));
    56                                 vxenum=i;
    57                         elseif (strcmpi(A(i).name(end),'y') | strncmpi(A(i).name,'vy',2));
    58                                 vyenum=i;
    59                         end
    60                 end
    61         end
     29%Interpolation
     30if strcmpi(Names.interp,'grid'),
     31        md.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,default_value);
     32        md.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,default_value);
    6233else
    63         error(['plugvelocities error message: file ' filename  ' not supported yet (it should hold 4 variables x,y,vx and vy (for grids) OR 5 variables  x,y,index,vx and vy (for mesh))']);
     34        md.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,default_value);
     35        md.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,default_value);
    6436end
    6537
    66 %assum that we have found at least vxenum and vyenum
    67 if ( isnan(vxenum) | isnan(vyenum))
    68         error(['plugvelocities error message: file ' filename  ' not supported yet (the velocities should be named vx and vy)']);
    69 end
    70 
    71 %find index
    72 if (~isgrid & insnan(indexenum)),
    73         for i=1:5
    74                 lengthi=min(A(i).size);
    75                 if (lengthi==3),
    76                         indexenum=i;
    77                 end
    78         end
    79         if isnan(indexenum),
    80                 error(['plugvelocities error message: file ' filename  ' not supported yet (index not found)']);
    81         end
    82 end
    83 
    84 %find x y
    85 if (isnan(xenum) | isnan(yenum))
    86 
    87         %check the size
    88         if A(vxenum).size(1)==A(vxenum).size(2),
    89                 error(['plugvelocities error message: file ' filename  ' not supported (velocities is a square matrix, save x and y with another name)']);
    90         end
    91         if ~(A(vxenum).size(1)==A(vyenum).size(1) & A(vxenum).size(2)==A(vyenum).size(2)),
    92                 error(['plugvelocities error message: file ' filename  ' not supported (vx and vy matrices do not have the same size)']);
    93         end
    94 
    95         %find xenum and yenum
    96         for i=1:4
    97                 lengthi=max(A(i).size);
    98                 if ((i~=vxenum) & (lengthi==A(vxenum).size(1) | lengthi==A(vxenum).size(1)+1)),
    99                         yenum=i;
    100                 elseif ((i~=vxenum) & (lengthi==A(vxenum).size(2) | lengthi==A(vxenum).size(2)+1)),
    101                         xenum=i;
    102                 end
    103         end
    104 
    105         %last check
    106         if (isnan(xenum) | isnan(yenum))
    107                 error(['plugdata error message: file ' filename  ' not supported yet']);
    108         end
    109 end
    110 
    111 %create names:
    112 xname=A(xenum).name;
    113 yname=A(yenum).name;
    114 vxname=A(vxenum).name;
    115 vyname=A(vyenum).name;
    116 if ~isgrid, indexname=A(indexenum).name; end
    117 
    118 %load data
    119 B=load(filename);
    120 
    121 %get x y vx and vy
    122 eval(['x=B.' xname ';']);
    123 eval(['y=B.' yname ';']);
    124 eval(['vx=B.' vxname ';']);
    125 eval(['vy=B.' vyname ';']);
    126 if ~isgrid, eval(['index=B.' indexname ';']); end
    127 
    128 %interpolate
    129 if isgrid,
    130         md.vx_obs=InterpFromGridToMesh(x,y,vx,md.x,md.y,default_value);
    131         md.vy_obs=InterpFromGridToMesh(x,y,vy,md.x,md.y,default_value);
    132 else
    133         md.vx_obs=InterpFromMeshToMesh2d(index,x,y,vx,md.x,md.y,default_value);
    134         md.vy_obs=InterpFromMeshToMesh2d(index,x,y,vy,md.x,md.y,default_value);
    135 end
     38md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    13639md.vx=md.vx_obs;
    13740md.vy=md.vy_obs;
    138 md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    13941md.vel=md.vel_obs;
Note: See TracChangeset for help on using the changeset viewer.