Changeset 515 for trunk/src/calc_field.m


Ignore:
Timestamp:
Aug 15, 2012, 11:36:12 PM (9 years ago)
Author:
sommeria
Message:

improvement of calc-field and combination of two fields

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/calc_field.m

    r501 r515  
    3838
    3939%% list of field options implemented
    40 FieldOptions={'velocity';...%image correlation corresponding to a vel vector
    41     'ima_cor';...%image correlation corresponding to a vel vector
    42     'norm_vel';...%norm of the velocity
    43     'vort';...%vorticity
    44     'div';...%divergence
    45     'strain';...%rate of strain
    46     'u';... %u velocity component
    47     'v';... %v velocity component
     40FieldOptions={'vec(U,V)';...%image correlation corresponding to a vel vector
     41    'C';...%image correlation corresponding to a vel vector
     42    'norm(U,V)';...%norm of the velocity
     43    'curl(U,V)';...%vorticity
     44    'div(U,V)';...%divergence
     45    'strain(U,V)';...%rate of strain
     46    'U';... %u velocity component
     47    'V';... %v velocity component
    4848    'w';... %w velocity component
    4949    'w_normal';... %w velocity component normal to the plane
     
    7777end
    7878FieldList=FieldList(check_calc==1);
    79 % if isempty(FieldList)
    80 %     DataOut=DataIn;
    81 %     return
    82 % end
    8379if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
    8480    nbcoord=3;
     
    9187units_cell={};
    9288
     89%% reproduce global attributes
     90DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;
     91for ilist=1:numel(DataOut.ListGlobalAttribute)
     92    DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
     93end
     94
    9395%% interpolation with new civ data
    94 if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
    95     %reproduce global attributes
    96     DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;
    97     for ilist=1:numel(DataOut.ListGlobalAttribute)
    98         DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
    99     end
    100 
     96[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn);
     97icell_tps=[];
     98for icell=1:numel(CellVarIndex)
     99    VarType=VarTypeCell{icell};
     100    if NbDimVec(icell)>=2 && ~isempty(VarType.subrange_tps)     
     101        icell_tps=[icell_tps icell];
     102    end
     103end
     104
     105%if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
     106if ~isempty(icell_tps)
    101107    %create a default grid if needed
    102108    if  ~exist('Coord_interp','var')
    103             XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
    104             YMax=max(max(DataIn.SubRange(2,:,:)));
    105             XMin=min(min(DataIn.SubRange(1,:,:)));
    106             YMin=min(min(DataIn.SubRange(2,:,:)));
     109        for ind=1:numel(icell_tps)
     110            SubRange=DataIn.(DataIn.ListVarName{VarType{icell_tps(ind)}.subrange_tps});
     111            XMax(ind)=max(max(SubRange(1,:,:)));% extrema of the coordinates
     112            YMax(ind)=max(max(SubRange(2,:,:)));
     113            XMin(ind)=min(min(SubRange(1,:,:)));
     114            YMin(ind)=min(min(SubRange(2,:,:)));
     115        end
     116        XMax=max(XMax);
     117        YMax=max(YMax);
     118        XMin=min(XMin);
     119        YMin=min(YMin);
    107120        if ~isfield(DataIn,'Mesh')
    108121            DataIn.Mesh=sqrt(2*(XMax-XMin)*(YMax-YMin)/numel(DataIn.Coord_tps));
     
    119132        coord_x=XMin:DataIn.Mesh:XMax;% increase the recommanded mesh to  visualisation
    120133        coord_y=YMin:DataIn.Mesh:YMax;
    121 %         npx=length(coord_x);
    122 %         npy=length(coord_y);
    123134        DataOut.coord_x=[XMin XMax];
    124135        DataOut.coord_y=[YMin YMax];
    125136        [XI,YI]=meshgrid(coord_x,coord_y);
    126 %         XI=reshape(XI,[],1);
    127 %         YI=reshape(YI,[],1);
    128137        Coord_interp=cat(3,XI,YI);%[XI YI];
    129138    end
     
    131140    npy=size(Coord_interp,1);
    132141    Coord_interp=reshape(Coord_interp,npx*npy,size(Coord_interp,3));
    133 %         npy=length(coord_y);
     142
    134143    %initialise output
    135144    nb_sites=size(Coord_interp,1);
     
    152161   
    153162    % interpolate data in each subdomain
    154     for isub=1:NbSubDomain
    155         nbvec_sub=DataIn.NbSites(isub);
    156         check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');
    157         ind_sel=find(sum(check_range,2)==nb_coord);
    158         %rho smoothing parameter
    159         %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
    160         %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
    161         nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
    162         if check_grid
    163             EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
    164         end
    165         if check_der
    166             [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'
    167         end
    168         ListVar={};
    169         for ilist=1:length(FieldList)
    170             var_count=numel(ListVar);
    171             switch FieldList{ilist}
    172                 case 'velocity'
    173                     ListVar=[ListVar {'U', 'V'}];
    174                     VarAttribute{var_count+1}.Role='vector_x';
    175                     VarAttribute{var_count+2}.Role='vector_y';
    176                     DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    177                     DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    178                 case 'u'
    179                     ListVar=[ListVar {'u'}];
    180                     VarAttribute{var_count+1}.Role='scalar';
    181                     DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    182                 case 'v'
    183                     ListVar=[ListVar {'v'}];
    184                     VarAttribute{var_count+1}.Role='scalar';
    185                     DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    186                 case 'norm_vel'
    187                     ListVar=[ListVar {'norm_vel'}];
    188                     VarAttribute{var_count+1}.Role='scalar';
    189                     U=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    190                     V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    191                     DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
    192                 case 'vort'
    193                     ListVar=[ListVar {'vort'}];
    194                     VarAttribute{var_count+1}.Role='scalar';
    195                     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);
    196                 case 'div'
    197                     ListVar=[ListVar {'div'}];
    198                     VarAttribute{var_count+1}.Role='scalar';
    199                     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);
    200                 case 'strain'
    201                     ListVar=[ListVar {'strain'}];
    202                     VarAttribute{var_count+1}.Role='scalar';
    203                     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);
     163    for icell=icell_tps
     164        for isub=1:NbSubDomain
     165            nbvec_sub=DataIn.(DataIn.ListVarName{VarType{icell}.nbsites_tps});
     166            SubRange=DataIn.(DataIn.ListVarName{VarType{icell}.subrange_tps});
     167            check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
     168            ind_sel=find(sum(check_range,2)==nb_coord);
     169            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     170            Coord_tps=DataIn.(DataIn.ListVarName{VarType{icell}.coord_tps});
     171            if check_grid
     172                EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
    204173            end
    205         end
    206     end
    207     DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
    208     nbval(nbval==0)=1;% to avoid division by zero for averaging
    209     if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
    210         DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
    211         DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
    212         DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
    213     end
    214     DataOut.ListVarName=[DataOut.ListVarName ListVar];
    215     for ifield=1:numel(ListVar)
    216         VarDimName{ifield}={'coord_y','coord_x'};
    217         DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
    218         DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
    219     end
    220     DataOut.FF=reshape(DataOut.FF,npy,npx);
    221     DataOut.VarDimName=[DataOut.VarDimName VarDimName];     
    222     DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
     174            if check_der
     175                [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
     176            end
     177            ListVar={};
     178            U_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_x_tps});
     179            V_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_y_tps});
     180            for ilist=1:length(FieldList)
     181                var_count=numel(ListVar);
     182                switch FieldList{ilist}
     183                    case 'velocity'
     184                        ListVar=[ListVar {'U', 'V'}];
     185                        VarAttribute{var_count+1}.Role='vector_x';
     186                        VarAttribute{var_count+2}.Role='vector_y';
     187                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
     188                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
     189                    case 'u'
     190                        ListVar=[ListVar {'u'}];
     191                        VarAttribute{var_count+1}.Role='scalar';
     192                        DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
     193                    case 'v'
     194                        ListVar=[ListVar {'v'}];
     195                        VarAttribute{var_count+1}.Role='scalar';
     196                        DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
     197                    case 'norm_vel'
     198                        ListVar=[ListVar {'norm_vel'}];
     199                        VarAttribute{var_count+1}.Role='scalar';
     200                        U=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
     201                        V=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
     202                        DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
     203                    case 'vort'
     204                        ListVar=[ListVar {'vort'}];
     205                        VarAttribute{var_count+1}.Role='scalar';
     206                        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);
     207                    case 'div'
     208                        ListVar=[ListVar {'div'}];
     209                        VarAttribute{var_count+1}.Role='scalar';
     210                        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);
     211                    case 'strain'
     212                        ListVar=[ListVar {'strain'}];
     213                        VarAttribute{var_count+1}.Role='scalar';
     214                        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);
     215                end
     216            end
     217        end
     218        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
     219        nbval(nbval==0)=1;% to avoid division by zero for averaging
     220        if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
     221            DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
     222            DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
     223            DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
     224        end
     225        DataOut.ListVarName=[DataOut.ListVarName ListVar];
     226        for ifield=1:numel(ListVar)
     227            VarDimName{ifield}={'coord_y','coord_x'};
     228            DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
     229            DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
     230        end
     231        DataOut.FF=reshape(DataOut.FF,npy,npx);
     232        DataOut.VarDimName=[DataOut.VarDimName VarDimName];
     233        DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
     234    end
    223235else
    224236
Note: See TracChangeset for help on using the changeset viewer.