Changeset 2574


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

better interpolation routine call

Location:
issm/trunk/src/m
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/classes/public/mesh/meshyams.m

    r2573 r2574  
    110110end
    111111if strcmpi(Names.interp,'grid'),
    112                 md.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
    113                 md.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
    114         else
    115                 md.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
    116                 md.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
    117         end
     112        md.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
     113        md.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
     114else
     115        md.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
     116        md.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
     117end
    118118md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);

  • 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;

  • issm/trunk/src/m/utils/Interp/InterpFromFile.m

    r2296 r2574  
    3535end
    3636
    37 %Get variables
    38 A=whos('-file',filename);
    39 
    40 %check the number of variables
    41 if (length(A)~=4 & length(A)~=3)
    42         error(['plugdata error message: file ' filename  ' not supported yet (it should hold 3 or 4 variables)']);
     37%load file
     38Names=FieldFindVarNames(filename);
     39Data=load(filename);
     40if strcmpi(Names.interp,'grid'),
     41        data_out=InterpFromGridToMesh(Data.(Names.xname),Data.(Names.yname),Data.(Names.dataname),md.x,md.y,default_value);
     42else
     43        data_out=InterpFromMeshToMesh2d(Data.(Names.indexname),Data.(Names.xname),Data.(Names.yname),Data.(Names.dataname),md.x,md.y,default_value);
    4344end
    44 
    45 %1: find data using their names
    46 xenum=NaN; yenum=NaN; indexenum=NaN; dataenum=NaN;
    47 for i=1:length(A),
    48 
    49         %x -> name begins with "x" and is a vector
    50         if (strncmpi(A(i).name,'x',1) & min(A(i).size)==1),
    51                 xenum=i;
    52 
    53         %y -> name begins with "y" and is a vector
    54         elseif (strncmpi(A(i).name,'y',1) & min(A(i).size)==1),
    55                 yenum=i;
    56 
    57         %index -> name begins with "index" or "elements" and 3 columns
    58         elseif length(A)==4
    59                 if ((strncmpi(A(i).name,'index',5) | strncmpi(A(i).name,'elements',8)) & A(i).size(2)==3),
    60                         indexenum=i;
    61                 end
    62 
    63         %data -> name begins with "field" and is a vector or a matrix
    64         %elseif (strcmpi(A(i).name,field) & max(A(i).size)>1),
    65         %       dataenum=i;
    66 
    67         end
    68 end
    69 
    70 %2: if only one item is missing, find it by elimination
    71 if length(A)==4
    72         pos=find(isnan([xenum yenum indexenum dataenum]));
    73         if length(pos)==1,
    74                 list=[xenum yenum indexenum dataenum]; list(pos)=[];
    75                 if pos==1,
    76                         xenum=setdiff(1:4,list);
    77                 elseif pos==2,
    78                         yenum=setdiff(1:4,list);
    79                 elseif pos==3,
    80                         indexenum=setdiff(1:4,list);
    81                 elseif pos==4,
    82                         dataenum=setdiff(1:4,list);
    83                 end
    84         end
    85 else
    86         pos=find(isnan([xenum yenum dataenum]));
    87         if length(pos)==1,
    88                 list=[xenum yenum indexenum dataenum]; list(pos)=[];
    89                 if pos==1,
    90                         xenum=setdiff(1:3,list);
    91                 elseif pos==2,
    92                         yenum=setdiff(1:3,list);
    93                 elseif pos==3,
    94                         dataenum=setdiff(1:3,list);
    95                 end
    96         end
    97 end
    98 
    99 %3: if one or several variables are missing, use sizes instead of names (works only if data is a matrix...)
    100 if (isnan(xenum) | isnan(yenum) |  isnan(dataenum) | (length(A)==4 & isnan(indexenum)))
    101 
    102         if (length(A)==4 & isnan(indexenum)),
    103                 for i=1:4
    104                         if A(i).size(2)==3,
    105                                 indexenum=i;
    106                         end
    107                 end
    108                 if isnan(indexenum)
    109                         error(['InterpFromFile error message: file ' filename  ' not supported yet (index not found)']);
    110                 end
    111         end
    112 
    113         if isnan(dataenum),
    114                 maxsize1=1;
    115                 maxsize2=1;
    116 
    117                 %find dataenum
    118                 for i=1:length(A)
    119                         sizei=A(i).size;
    120                         if (sizei(1)>=maxsize1-1 & sizei(2)>=maxsize2-1),
    121                                 dataenum=i
    122                                 maxsize1=sizei(1);
    123                                 maxsize2=sizei(2);
    124                         end
    125                 end
    126                 if isnan(dataenum)
    127                         error(['InterpFromFile error message: file ' filename  ' not supported yet (data not found)']);
    128                 end
    129                 if (A(dataenum).size(1)==A(dataenum).size(2) & isnan(xenum) & isnan(yenum)),
    130                         error(['InterpFromFile error message: file ' filename  ' not supported (the data is a square matrix, save x and y with another name)']);
    131                 end
    132         end
    133 
    134         %find xenum and yenum
    135         if isnan(xenum) | isnan(yenum)
    136                 for i=1:length(A),
    137                         lengthi=max(A(i).size);
    138                         if isnan(yenum) & ((i~=dataenum) & (lengthi==A(dataenum).size(1) | lengthi==A(dataenum).size(1)+1)),
    139                                 yenum=i;
    140                         elseif isnan(xenum) & ((i~=dataenum) & (lengthi==A(dataenum).size(2) | lengthi==A(dataenum).size(2)+1)),
    141                                 xenum=i;
    142                         end
    143                 end
    144         end
    145 
    146         %last chance: if only one item is missing, find it by elimination
    147         if length(A)==4
    148                 pos=find(isnan([xenum yenum indexenum dataenum]));
    149                 if length(pos)==1,
    150                         list=[xenum yenum indexenum dataenum]; list(pos)=[];
    151                         if pos==1,
    152                                 xenum=setdiff(1:4,list);
    153                         elseif pos==2,
    154                                 yenum=setdiff(1:4,list);
    155                         elseif pos==3,
    156                                 indexenum=setdiff(1:4,list);
    157                         elseif pos==4,
    158                                 dataenum=setdiff(1:4,list);
    159                         end
    160                 end
    161         else
    162                 pos=find(isnan([xenum yenum dataenum]));
    163                 if length(pos)==1,
    164                         list=[xenum yenum indexenum dataenum]; list(pos)=[];
    165                         if pos==1,
    166                                 xenum=setdiff(1:3,list);
    167                         elseif pos==2,
    168                                 yenum=setdiff(1:3,list);
    169                         elseif pos==3,
    170                                 dataenum=setdiff(1:3,list);
    171                         end
    172                 end
    173         end
    174 
    175         %last check
    176         if (isnan(xenum) | isnan(yenum) | isnan(dataenum) | (length(A)==4 & isnan(indexenum)))
    177                 error(['InterpFromFile error message: file ' filename  ' not supported yet (not recognzed variable names or sizes)']);
    178         end
    179 end
    180 
    181 %OK, interpolate
    182 if length(A)==4,
    183 
    184         %create names:
    185         xname=A(xenum).name;
    186         yname=A(yenum).name;
    187         indexname=A(indexenum).name;
    188         dataname=A(dataenum).name;
    189 
    190         %load data
    191         B=load(filename);
    192 
    193         %get x y and data
    194         eval(['x_data=B.' xname ';'])
    195         eval(['y_data=B.' yname ';'])
    196         eval(['index_data=B.' indexname ';'])
    197         eval(['data=B.' dataname ';'])
    198 
    199         %interpolate
    200         data_out=InterpFromMeshToMesh2d(index_data,x_data(:),y_data(:),data(:),x(:),y(:),default_value);
    201 
    202 else
    203 
    204         %create names:
    205         xname=A(xenum).name;
    206         yname=A(yenum).name;
    207         dataname=A(dataenum).name;
    208 
    209         %load data
    210         B=load(filename);
    211 
    212         %get x y and data
    213         eval(['x_data=B.' xname ';'])
    214         eval(['y_data=B.' yname ';'])
    215         eval(['data=B.' dataname ';'])
    216 
    217         %interpolate
    218         data_out=InterpFromGridToMesh(x_data(:),y_data(:),data,x(:),y(:),default_value);
    219 end

  • issm/trunk/src/m/utils/Interp/VelFindVarNames.m

    r2573 r2574  
    11function Names=VelFindVarNames(filename)
    22%VELFINDVARNAMES - find names of variables in a velocity data set file
     3%
     4%   This routines looks at the variables contained in a file and finds out
     5%   the names of the variables that are needed for an interpolation (x,y,vx,vy)
     6%   or (index,x,y,vx,vy)
    37%
    48%   Usage:
Note: See TracChangeset for help on using the changeset viewer.