Changeset 404 for trunk/src/calc_field.m
 Timestamp:
 Apr 30, 2012, 6:46:45 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/calc_field.m
r399 r404 1 2 1 3 %'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them 2 4 % … … 35 37 % the scalar is calculated from other fields, as explicited below 36 38 37 list_field={'velocity';...%image correlation corresponding to a vel vector 39 %% list of field options implemented 40 FieldOptions={'velocity';...%image correlation corresponding to a vel vector 38 41 'ima_cor';...%image correlation corresponding to a vel vector 39 42 'norm_vel';...%norm of the velocity … … 48 51 errormsg=[]; %default error message 49 52 if ~exist('FieldList','var') 50 DataOut=list_field;% gives the list of possible fields in the absence of input 53 DataOut=FieldOptions;% gives the list of possible field inputs in the absence of input 54 return 55 end 56 57 %% check input 58 if ~exist('DataIn','var') 59 DataIn=[]; 60 end 61 if ischar(FieldList) 62 FieldList={FieldList};%convert a string input to a cell with one string element 63 end 64 check_grid=0; 65 check_der=0; 66 check_calc=ones(size(FieldList)); 67 for ilist=1:length(FieldList) 68 switch FieldList{ilist} 69 case {'u','v'} 70 check_grid=1;% needs a regular grid 71 case{'vort','div','strain'}% needs spatial derivatives spatial derivatives 72 check_der=1; 73 case {'velocity','norm_vel'}; 74 otherwise 75 check_calc(ilist)=0; 76 end 77 end 78 FieldList=FieldList(check_calc==1); 79 if isempty(FieldList) 80 DataOut=DataIn; 81 return 82 end 83 if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X)) 84 nbcoord=3; 51 85 else 52 if ~exist('DataIn','var') 53 DataIn=[]; 54 end 55 if ischar(FieldList) 56 FieldList={FieldList};%convert a string input to a cell with one string element 57 end 58 if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X)) 59 nbcoord=3; 60 else 61 nbcoord=2; 62 end 63 ListVarName={}; 64 ValueList={}; 65 RoleList={}; 66 units_cell={}; 67 68 %% interpolation with new civ data 69 if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps')&& exist('Coord_interp','var') 70 DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute 71 for ilist=1:numel(DataOut.ListGlobalAttribute) 72 DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist}); 73 end 74 DataOut.ListVarName={'coord_y','coord_x','FF'}; 75 DataOut.VarDimName{1}='coord_y'; 76 DataOut.VarDimName{2}='coord_x'; 77 XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates 78 YMax=max(max(DataIn.SubRange(2,:,:))); 79 XMin=min(min(DataIn.SubRange(1,:,:))); 80 YMin=min(min(DataIn.SubRange(2,:,:))); 81 check_der=0; 82 check_val=0; 83 nb_sites=size(Coord_interp,1); 84 nb_coord=size(Coord_interp,2); 86 nbcoord=2; 87 end 88 ListVarName={}; 89 ValueList={}; 90 RoleList={}; 91 units_cell={}; 92 93 %% interpolation with new civ data 94 if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var')  check_grid check_der) 95 %create a default grid if needed 96 if ~exist('Coord_interp','var') 97 coord_x=DataIn.XMin:DataIn.Mesh:DataIn.XMax; 98 coord_y=DataIn.XMin:DataIn.Mesh:DataIn.YMax; 99 [XI,YI]=meshgrid(coord_x,coord_y); 100 XI=reshape(XI,[],1); 101 YI=reshape(YI,[],1); 102 Coord_interp=[XI YI]; 103 end 104 DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute 105 for ilist=1:numel(DataOut.ListGlobalAttribute) 106 DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist}); 107 end 108 DataOut.ListVarName={'coord_y','coord_x','FF'}; 109 DataOut.VarDimName{1}='coord_y'; 110 DataOut.VarDimName{2}='coord_x'; 111 XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates 112 YMax=max(max(DataIn.SubRange(2,:,:))); 113 XMin=min(min(DataIn.SubRange(1,:,:))); 114 YMin=min(min(DataIn.SubRange(2,:,:))); 115 % check_der=0; 116 % check_val=0; 117 nb_sites=size(Coord_interp,1); 118 nb_coord=size(Coord_interp,2); 119 %initialise output 120 for ilist=1:length(FieldList) 121 switch FieldList{ilist} 122 case 'velocity' 123 DataOut.U=zeros(nb_sites,1); 124 DataOut.V=zeros(nb_sites,1); 125 otherwise 126 % case{'vort','div','strain'}% case of spatial derivatives 127 DataOut.(FieldList{ilist})=zeros(nb_sites,1); 128 % otherwise % case of a scalar 129 % check_val=1; 130 % DataOut.(FieldList{ilist})=zeros(size(Coord_interp,1)); 131 end 132 end 133 nbval=zeros(nb_sites,1); 134 NbSubDomain=size(DataIn.SubRange,3); 135 %DataIn.Coord_tps=DataIn.Coord_tps(1:end3,:,:);% suppress the 3 zeros used to fit with the dimensions of variables 136 for isub=1:NbSubDomain 137 nbvec_sub=DataIn.NbSites(isub); 138 check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)'); 139 ind_sel=find(sum(check_range,2)==nb_coord); 140 %rho smoothing parameter 141 % epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites 142 % ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs 143 nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap) 144 if check_val 145 EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources' 146 end 147 if check_der 148 [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' 149 end 85 150 for ilist=1:length(FieldList) 86 151 switch FieldList{ilist} 87 152 case 'velocity' 88 check_val=1; 89 DataOut.U=zeros(nb_sites,1); 90 DataOut.V=zeros(nb_sites,1); 91 case{'vort','div','strain'}% case of spatial derivatives 92 check_der=1; 93 DataOut.(FieldList{ilist})=zeros(nb_sites,1); 94 % otherwise % case of a scalar 95 % check_val=1; 96 % DataOut.(FieldList{ilist})=zeros(size(Coord_interp,1)); 153 ListFields={'U', 'V'}; 154 VarAttributes{1}.Role='vector_x'; 155 VarAttributes{2}.Role='vector_y'; 156 DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); 157 DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); 158 case 'u' 159 ListFields={'U'}; 160 VarAttributes{1}.Role='scalar'; 161 DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); 162 case 'v' 163 ListFields={'V'}; 164 VarAttributes{1}.Role='scalar'; 165 DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); 166 case 'norm_vel' 167 ListFields={'norm_vel'}; 168 VarAttributes{1}.Role='scalar'; 169 V=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); 170 V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); 171 DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V); 172 case 'vort' 173 ListFields={'vort'}; 174 VarAttributes{1}.Role='scalar'; 175 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); 176 case 'div' 177 ListFields={'div'}; 178 VarAttributes{1}.Role='scalar'; 179 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); 180 case 'strain' 181 ListFields={'strain'}; 182 VarAttributes{1}.Role='scalar'; 183 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); 97 184 end 98 185 end 99 nbval=zeros(nb_sites,1); 100 NbSubDomain=size(DataIn.SubRange,3); 101 %DataIn.Coord_tps=DataIn.Coord_tps(1:end3,:,:);% suppress the 3 zeros used to fit with the dimensions of variables 102 for isub=1:NbSubDomain 103 nbvec_sub=DataIn.NbSites(isub); 104 check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)'); 105 ind_sel=find(sum(check_range,2)==nb_coord); 106 %rho smoothing parameter 107 % epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites 108 % ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs 109 nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap) 110 if check_val 111 EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources' 112 end 113 if check_der 114 [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' 115 end 116 for ilist=1:length(FieldList) 117 switch FieldList{ilist} 118 case 'velocity' 119 ListFields={'U', 'V'}; 120 VarAttributes{1}.Role='vector_x'; 121 VarAttributes{2}.Role='vector_y'; 122 DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); 123 DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); 124 case 'u' 125 ListFields={'U'}; 126 VarAttributes{1}.Role='scalar'; 127 DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); 128 case 'v' 129 ListFields={'V'}; 130 VarAttributes{1}.Role='scalar'; 131 DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); 132 case 'vort' 133 ListFields={'vort'}; 134 VarAttributes{1}.Role='scalar'; 135 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); 136 case 'div' 137 ListFields={'div'}; 138 VarAttributes{1}.Role='scalar'; 139 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); 140 case 'strain' 141 ListFields={'strain'}; 142 VarAttributes{1}.Role='scalar'; 143 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); 144 end 145 end 146 DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang 147 % DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI)); 148 nbval(nbval==0)=1; 149 % switch FieldList{1} 150 % case {'velocity','u','v'} 151 % DataOut.U=reshape(DataOut.U./nbval,numel(yI),numel(xI)); 152 % DataOut.V=reshape(DataOut.V./nbval,numel(yI),numel(xI)); 153 % case 'vort' 154 % DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI)); 155 % case 'div' 156 % DataOut.div=reshape(DataOut.div,numel(yI),numel(xI)); 157 % case 'strain' 158 % DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI)); 159 % end 160 DataOut.ListVarName=[DataOut.ListVarName ListFields]; 161 for ilist=3:numel(DataOut.ListVarName) 162 DataOut.VarDimName{ilist}={'coord_y','coord_x'}; 163 end 164 DataOut.VarAttribute={[],[]}; 165 DataOut.VarAttribute{3}.Role='errorflag'; 166 DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes]; 167 end 168 else 169 170 %% civx data 171 DataOut=DataIn; 172 for ilist=1:length(FieldList) 173 if ~isempty(FieldList{ilist}) 174 [VarName,Value,Role,units]=feval(FieldList{ilist},DataIn);%calculate field with appropriate function named FieldList{ilist} 175 ListVarName=[ListVarName VarName]; 176 ValueList=[ValueList Value]; 177 RoleList=[RoleList Role]; 178 units_cell=[units_cell units]; 179 end 180 end 181 %erase previous data (except coordinates) 182 for ivar=nbcoord+1:length(DataOut.ListVarName) 183 VarName=DataOut.ListVarName{ivar}; 184 DataOut=rmfield(DataOut,VarName); 185 end 186 DataOut.ListVarName=DataOut.ListVarName(1:nbcoord); 187 if isfield(DataOut,'VarDimName') 188 DataOut.VarDimName=DataOut.VarDimName(1:nbcoord); 189 else 190 errormsg='element .VarDimName missing in input data'; 191 return 192 end 193 DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord); 194 %append new data 195 DataOut.ListVarName=[DataOut.ListVarName ListVarName]; 196 for ivar=1:length(ListVarName) 197 DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1}; 198 DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar}; 199 DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar}; 200 DataOut.(ListVarName{ivar})=ValueList{ivar}; 201 end 202 end 203 end 186 DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang 187 nbval(nbval==0)=1; 188 DataOut.ListVarName=[DataOut.ListVarName ListFields]; 189 for ilist=3:numel(DataOut.ListVarName) 190 DataOut.VarDimName{ilist}={'coord_y','coord_x'}; 191 end 192 DataOut.VarAttribute={[],[]}; 193 DataOut.VarAttribute{3}.Role='errorflag'; 194 DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes]; 195 end 196 else 197 198 %% civx data 199 DataOut=DataIn; 200 for ilist=1:length(FieldList) 201 if ~isempty(FieldList{ilist}) 202 [VarName,Value,Role,units]=feval(FieldList{ilist},DataIn);%calculate field with appropriate function named FieldList{ilist} 203 ListVarName=[ListVarName VarName]; 204 ValueList=[ValueList Value]; 205 RoleList=[RoleList Role]; 206 units_cell=[units_cell units]; 207 end 208 end 209 %erase previous data (except coordinates) 210 for ivar=nbcoord+1:length(DataOut.ListVarName) 211 VarName=DataOut.ListVarName{ivar}; 212 DataOut=rmfield(DataOut,VarName); 213 end 214 DataOut.ListVarName=DataOut.ListVarName(1:nbcoord); 215 if isfield(DataOut,'VarDimName') 216 DataOut.VarDimName=DataOut.VarDimName(1:nbcoord); 217 else 218 errormsg='element .VarDimName missing in input data'; 219 return 220 end 221 DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord); 222 %append new data 223 DataOut.ListVarName=[DataOut.ListVarName ListVarName]; 224 for ivar=1:length(ListVarName) 225 DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1}; 226 DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar}; 227 DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar}; 228 DataOut.(ListVarName{ivar})=ValueList{ivar}; 229 end 230 end 231 204 232 205 233
Note: See TracChangeset
for help on using the changeset viewer.