Changeset 397 for trunk/src/calc_field.m


Ignore:
Timestamp:
Apr 26, 2012, 8:59:09 AM (12 years ago)
Author:
sommeria
Message:

civ_matlab and patch improved, changes in the management of interpolation (still in progress).
adapatation to movies (use of VideoReader?)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/calc_field.m

    r389 r397  
    2323% file, depending on the scalar
    2424
    25 function [DataOut,errormsg]=calc_field(FieldName,DataIn,VelType,XI,YI)
     25function [DataOut,errormsg]=calc_field(FieldName,DataIn,Coord_interp)
    2626
    2727%list of defined scalars to display in menus (in addition to 'ima_cor').
    … …  
    6666   
    6767    %% interpolation with new civ data
    68     if isfield(DataIn,'SubRange') &&(strcmp(VelType,'filter1')||strcmp(VelType,'filter2'))
     68    if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps')&& exist('Coord_interp','var')
     69        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
     70        for ilist=1:numel(DataOut.ListGlobalAttribute)
     71            DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
     72        end
     73        DataOut.ListVarName={'coord_y','coord_x','FF'};
     74        DataOut.VarDimName{1}='coord_y';
     75        DataOut.VarDimName{2}='coord_x';
    6976        XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
    7077        YMax=max(max(DataIn.SubRange(2,:,:)));
    7178        XMin=min(min(DataIn.SubRange(1,:,:)));
    7279        YMin=min(min(DataIn.SubRange(2,:,:)));
    73         if ~(exist('XI','var')&&exist('YI','var'))
    74             Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.Coord_tps(:,1,:)));%2D
    75             xI=XMin:Mesh:XMax;
    76             yI=YMin:Mesh:YMax;
    77             [XI,YI]=meshgrid(xI,yI);% interpolation grid
    78             XI=reshape(XI,[],1);
    79             YI=reshape(YI,[],1);
    80         end
    81         DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
    82         for ilist=1:numel(DataOut.ListGlobalAttribute)
    83             DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
    84         end
    85         DataOut.ListVarName={'coord_y','coord_x','FF'};
    86         DataOut.VarDimName{1}='coord_y';
    87         DataOut.VarDimName{2}='coord_x';
    88         DataOut.coord_y=[yI(1) yI(end)];
    89         DataOut.coord_x=[xI(1) xI(end)];
    90         switch FieldName{1}
    91             case {'velocity','u','v'}
    92                 DataOut.U=zeros(size(XI));
    93                 DataOut.V=zeros(size(XI));
    94             case 'vort'
    95                 DataOut.vort=zeros(size(XI));
    96             case 'div'
    97                 DataOut.div=zeros(size(XI));
    98             case 'strain'
    99                 DataOut.strain=zeros(size(XI));
    100         end
    101         nbval=zeros(size(XI));
     80        %         if ~exist('Coord_interp','var')
     81        %             Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.Coord_tps(:,1,:)));%2D
     82        %             xI=XMin:Mesh:XMax;
     83        %             yI=YMin:Mesh:YMax;
     84        %             [XI,YI]=meshgrid(xI,yI);% interpolation grid
     85        %             Coord_interp(:,1)=reshape(XI,[],1);
     86        %             Coord_interp(:,2)=reshape(YI,[],1);
     87        %             DataOut.coord_y=[yI(1) yI(end)];
     88        %             DataOut.coord_x=[xI(1) xI(end)];
     89        %      end
     90        check_der=0;
     91        check_val=0;
     92        for ilist=1:length(FieldName)
     93            switch FieldName{ilist}
     94                case 'velocity'
     95                    check_val=1;
     96                    DataOut.U=zeros(size(Coord_interp,1));
     97                    DataOut.V=zeros(size(Coord_interp,1));
     98                case{'vort','div','strain'}% case of spatial derivatives
     99                    check_der=1;
     100                    DataOut.(FieldName{ilist})=zeros(size(Coord_interp,1));
     101                otherwise % case of a scalar
     102                    check_val=1;
     103                    DataOut.(FieldName{ilist})=zeros(size(Coord_interp,1));
     104            end
     105        end
     106        nbval=zeros(size(Coord_interp,1));
    102107        NbSubDomain=size(DataIn.SubRange,3);
    103108        for isub=1:NbSubDomain
    … …  
    107112                nbvec_sub=numel(find(DataIn.Indices_tps(:,isub)));
    108113            end
    109             ind_sel=find(XI>=DataIn.SubRange(1,1,isub) & XI<=DataIn.SubRange(1,2,isub) & YI>=DataIn.SubRange(2,1,isub) & YI<=DataIn.SubRange(2,2,isub));
     114            ind_sel=find(Coord_interp>=DataIn.SubRange(:,1,isub) & Coord_interp<=DataIn.SubRange(:,2,isub));
    110115            %rho smoothing parameter
    111             epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites
    112             ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
     116            %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
     117            %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
    113118            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     119            if check_val
     120                EM = tps_eval(Coord_interp(ind_sel),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
     121            end
     122            if check_der
     123                [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
     124            end
     125            for ilist=1:length(FieldName)
     126                switch FieldName{ilist}
     127                    case 'velocity'
     128                        ListFields={'U', 'V'};
     129                        VarAttributes{1}.Role='vector_x';
     130                        VarAttributes{2}.Role='vector_y';
     131                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
     132                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     133                    case 'u'
     134                        ListFields={'U'};
     135                        VarAttributes{1}.Role='scalar';
     136                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
     137                    case 'v'
     138                        ListFields={'V'};
     139                        VarAttributes{1}.Role='scalar';
     140                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     141                    case 'vort'
     142                        ListFields={'vort'};
     143                        VarAttributes{1}.Role='scalar';
     144                        DataOut.vort(ind_sel)=DataOut.vort(ind_sel)+EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)-EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
     145                    case 'div'
     146                        ListFields={'div'};
     147                        VarAttributes{1}.Role='scalar';
     148                        DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
     149                    case 'strain'
     150                        ListFields={'strain'};
     151                        VarAttributes{1}.Role='scalar';
     152                        DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
     153                end
     154            end
     155            DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
     156            DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
     157            nbval(nbval==0)=1;
    114158            switch FieldName{1}
    115159                case {'velocity','u','v'}
    116                     EM = tps_eval(epoints,ctrs);%kernels for calculating the velocity from tps 'sources'
    117                 case{'vort','div','strain'}
    118                     [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%kernels for calculating the spatial derivatives from tps 'sources'
    119             end
    120             switch FieldName{1}
    121                 case 'velocity'
    122                     ListFields={'U', 'V'};
    123                     VarAttributes{1}.Role='vector_x';
    124                     VarAttributes{2}.Role='vector_y';
    125                     DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    126                     DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    127                 case 'u'
    128                     ListFields={'U'};
    129                     VarAttributes{1}.Role='scalar';
    130                     DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    131                 case 'v'
    132                     ListFields={'V'};
    133                     VarAttributes{1}.Role='scalar';
    134                     DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     160                    DataOut.U=reshape(DataOut.U./nbval,numel(yI),numel(xI));
     161                    DataOut.V=reshape(DataOut.V./nbval,numel(yI),numel(xI));
    135162                case 'vort'
    136                     ListFields={'vort'};
    137                     VarAttributes{1}.Role='scalar';
    138                     DataOut.vort(ind_sel)=DataOut.vort(ind_sel)+EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)-EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
     163                    DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
    139164                case 'div'
    140                     ListFields={'div'};
    141                     VarAttributes{1}.Role='scalar';
    142                     DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
     165                    DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
    143166                case 'strain'
    144                     ListFields={'strain'};
    145                     VarAttributes{1}.Role='scalar';
    146                     DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
    147             end
    148         end
    149         DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang     
    150         DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
    151         nbval(nbval==0)=1;
    152         switch FieldName{1}
    153               case {'velocity','u','v'}
    154         DataOut.U=reshape(DataOut.U./nbval,numel(yI),numel(xI));
    155         DataOut.V=reshape(DataOut.V./nbval,numel(yI),numel(xI));
    156             case 'vort'
    157         DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
    158             case 'div'
    159         DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
    160             case 'strain'
    161            DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI));
    162         end
    163         DataOut.ListVarName=[DataOut.ListVarName ListFields];
    164         for ilist=3:numel(DataOut.ListVarName)
    165             DataOut.VarDimName{ilist}={'coord_y','coord_x'};
    166         end
    167         DataOut.VarAttribute={[],[]};
    168         DataOut.VarAttribute{3}.Role='errorflag';
    169         DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
     167                    DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI));
     168            end
     169            DataOut.ListVarName=[DataOut.ListVarName ListFields];
     170            for ilist=3:numel(DataOut.ListVarName)
     171                DataOut.VarDimName{ilist}={'coord_y','coord_x'};
     172            end
     173            DataOut.VarAttribute={[],[]};
     174            DataOut.VarAttribute{3}.Role='errorflag';
     175            DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
     176        end
    170177    else
    171178       
    … …  
    200207            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
    201208            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
    202             eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
     209            DataOut.(ListVarName{ivar})=ValueList{ivar};
    203210        end
    204211    end
Note: See TracChangeset for help on using the changeset viewer.