Changeset 397 for trunk/src/calc_field.m
 Timestamp:
 Apr 26, 2012, 8:59:09 AM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/calc_field.m
r389 r397 23 23 % file, depending on the scalar 24 24 25 function [DataOut,errormsg]=calc_field(FieldName,DataIn, VelType,XI,YI)25 function [DataOut,errormsg]=calc_field(FieldName,DataIn,Coord_interp) 26 26 27 27 %list of defined scalars to display in menus (in addition to 'ima_cor'). … … 66 66 67 67 %% 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'; 69 76 XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates 70 77 YMax=max(max(DataIn.SubRange(2,:,:))); 71 78 XMin=min(min(DataIn.SubRange(1,:,:))); 72 79 YMin=min(min(DataIn.SubRange(2,:,:))); 73 if ~(exist('XI','var')&&exist('YI','var')) 74 Mesh=sqrt((YMaxYMin)*(XMaxXMin)/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((YMaxYMin)*(XMaxXMin)/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)); 102 107 NbSubDomain=size(DataIn.SubRange,3); 103 108 for isub=1:NbSubDomain … … 107 112 nbvec_sub=numel(find(DataIn.Indices_tps(:,isub))); 108 113 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)); 110 115 %rho smoothing parameter 111 epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites112 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs116 % epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites 117 % ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs 113 118 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; 114 158 switch FieldName{1} 115 159 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)); 135 162 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)); 139 164 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)); 143 166 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 170 177 else 171 178 … … 200 207 DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar}; 201 208 DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar}; 202 eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])209 DataOut.(ListVarName{ivar})=ValueList{ivar}; 203 210 end 204 211 end
Note: See TracChangeset
for help on using the changeset viewer.