Changeset 389 for trunk/src/calc_field.m


Ignore:
Timestamp:
Apr 8, 2012, 11:11:38 PM (12 years ago)
Author:
sommeria
Message:

several bugs corrected: object managing, tps filter...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/calc_field.m

    r383 r389  
    6565    units_cell={};
    6666   
    67     %% interpolation with new civ data
    68     %if isfield(DataIn,'X_SubRange')
     67    %% interpolation with new civ data
    6968    if isfield(DataIn,'SubRange') &&(strcmp(VelType,'filter1')||strcmp(VelType,'filter2'))
    70         XMax=max(max(DataIn.SubRange(1,:,:)));
     69        XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
    7170        YMax=max(max(DataIn.SubRange(2,:,:)));
    7271        XMin=min(min(DataIn.SubRange(1,:,:)));
    … …  
    7675            xI=XMin:Mesh:XMax;
    7776            yI=YMin:Mesh:YMax;
    78             [XI,YI]=meshgrid(xI,yI);
     77            [XI,YI]=meshgrid(xI,yI);% interpolation grid
    7978            XI=reshape(XI,[],1);
    8079            YI=reshape(YI,[],1);
    … …  
    8281        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
    8382        for ilist=1:numel(DataOut.ListGlobalAttribute)
    84             eval(['DataOut.' DataOut.ListGlobalAttribute{ilist} '=DataIn.' DataIn.ListGlobalAttribute{ilist} ';'])
     83            DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
    8584        end
    8685        DataOut.ListVarName={'coord_y','coord_x','FF'};
    … …  
    8988        DataOut.coord_y=[yI(1) yI(end)];
    9089        DataOut.coord_x=[xI(1) xI(end)];
    91         DataOut.U=zeros(size(XI));
    92         DataOut.V=zeros(size(XI));
    93         DataOut.vort=zeros(size(XI));       
    94         DataOut.div=zeros(size(XI));
     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
    95101        nbval=zeros(size(XI));
    96102        NbSubDomain=size(DataIn.SubRange,3);
    97103        for isub=1:NbSubDomain
    98104            if isfield(DataIn,'NbSites')
    99                 nbvec_sub=DataIn.NbSites(isub); 
     105                nbvec_sub=DataIn.NbSites(isub);
    100106            else
    101                 nbvec_sub=numel(find(DataIn.Indices_tps(:,isub))); 
     107                nbvec_sub=numel(find(DataIn.Indices_tps(:,isub)));
    102108            end
    103                 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));
    104                 %rho smoothing parameter
    105                 epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites
    106                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
    107                 nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
    108 %                 PM = [ones(size(epoints,1),1) epoints];
    109                 switch FieldName{1}
    110                     case {'velocity','u','v'}
    111                        % DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
    112                        % EM = tps(1,DM_eval);%values of thin plate
    113                         EM = tps_eval(epoints,ctrs);
    114                     case{'vort','div'}
    115                         [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
    116      
    117                 end
    118                 switch FieldName{1}
    119                     case 'velocity'
    120                         ListFields={'U', 'V'};
    121                         VarAttributes{1}.Role='vector_x';
    122                         VarAttributes{2}.Role='vector_y';
    123                         DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    124                         DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    125                     case 'u'
    126                         ListFields={'U'};
    127                         VarAttributes{1}.Role='scalar';
    128                         DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    129                     case 'v'
    130                         ListFields={'V'};
    131                         VarAttributes{1}.Role='scalar';
    132                         DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    133                     case 'vort'
    134                         ListFields={'vort'};
    135                         VarAttributes{1}.Role='scalar';
    136                         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);
    137                     case 'div'
    138                         ListFields={'div'};
    139                         VarAttributes{1}.Role='scalar';
    140                         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);
    141                 end
    142         end
    143         DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang 
    144        
     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));
     110            %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
     113            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     114            switch FieldName{1}
     115                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);
     135                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);
     139                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);
     143                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     
    145150        DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
    146151        nbval(nbval==0)=1;
    147         DataOut.U=DataOut.U ./nbval;
    148         DataOut.V=DataOut.V ./nbval;
    149         DataOut.U=reshape(DataOut.U,numel(yI),numel(xI));
    150         DataOut.V=reshape(DataOut.V,numel(yI),numel(xI));
     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'
    151157        DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
     158            case 'div'
    152159        DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
     160            case 'strain'
     161           DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI));
     162        end
    153163        DataOut.ListVarName=[DataOut.ListVarName ListFields];
    154164        for ilist=3:numel(DataOut.ListVarName)
    … …  
    157167        DataOut.VarAttribute={[],[]};
    158168        DataOut.VarAttribute{3}.Role='errorflag';
    159         DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];   
    160     else   
     169        DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
     170    else
    161171       
    162     %% civx data   
     172        %% civx data
    163173        DataOut=DataIn;
    164174        for ilist=1:length(FieldName)
Note: See TracChangeset for help on using the changeset viewer.