Changeset 1154 for trunk/src/filter_tps.m
 Timestamp:
 Jul 7, 2024, 11:22:00 PM (3 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/filter_tps.m
r1152 r1154 6 6 % SubRange(NbCoord,2,NbSubdomain): range (min, max) of the coordinates x and y respectively, for each subdomain 7 7 % NbCentre(NbSubdomain): number of source points for each subdomain 8 % FF: false flags preserved from the input, or equal to 20for vectors excluded by the difference with the smoothed field8 % FF: false flags preserved from the input, or equal to true for vectors excluded by the difference with the smoothed field 9 9 % U_smooth, V_smooth: filtered velocity components at the positions of the initial data 10 10 % Coord_tps(NbCentre,NbCoord,NbSubdomain): positions of the tps centres … … 61 61 62 62 %% smoothing parameter: CHANGED 03 May 2024 TO GET RESULTS INDEPENDENT OF SUBDOMAINSIZE 63 %smoothing=Siz(1)*Siz(2)*FieldSmooth/1000%o ptimum smoothing increase as the area of the subdomain (division by 1000 to reach good values with the default GUI input)63 %smoothing=Siz(1)*Siz(2)*FieldSmooth/1000%old calculation before 03 May < r1129 64 64 NbVecSub=NbVec/NbSubDomain;% refined estimation of the nbre of vectors per subdomain 65 65 smoothing=sqrt(Siz(1)*Siz(2)/NbVecSub)*FieldSmooth;%optimum smoothing increase as the typical mesh size =sqrt(SizX*SizY/NbVecSub)^1/2 66 Threshold=Threshold*Threshold;% take the square of the threshold to work with the modulus squared (not done before r1154) 66 67 67 68 %% default output … … 75 76 V_smooth=zeros(NbVec,1);% smoothed velocity V at the initial positions 76 77 W_smooth=[];%default (2 component case) 77 FF= zeros(NbVec,1);78 FF=false(NbVec,1);%false flag=0 (false) by default 78 79 nb_select=zeros(NbVec,1); 79 80 check_empty=zeros(1,NbSubDomain); … … 88 89 while numel(ind_sel)>numel(ind_sel_previous) 89 90 ind_sel_previous=ind_sel;% record the set of selected vector indices for next iteration 90 ind_sel= find(FF==0 & Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));% indices of vectors in the subdomain #isub 91 ind_sel= find(~FF & Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));% indices of vectors in the subdomain #isub 92 %ind_sel= find( Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));% indices of vectors in the subdomain #isub 91 93 % if no vector in the subdomain #isub, skip the subdomain 92 94 if isempty(ind_sel) … … 97 99 SubRange(:,1,isub)=SubRange(:,1,isub)Siz'/4; 98 100 SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4; 99 % subdomain includes enough vectors, perform tps interpolation101 % if subdomain includes enough vectors, perform tps interpolation 100 102 else 101 103 [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel,:),U(ind_sel),smoothing); … … 106 108 ind_ind_sel=1:numel(ind_sel);%default 107 109 if exist('Threshold','var')&&~isempty(Threshold) 108 FF(ind_sel)= 4*(NormDiff>Threshold);%put FF value to 4to identify the criterium of elimmination109 ind_ind_sel=find( FF(ind_sel)==0); % select the indices of remaining vectors in the subset of ind_sel vectors110 FF(ind_sel)=(NormDiff>Threshold);%put FF value to 1 to identify the criterium of elimmination 111 ind_ind_sel=find(~FF(ind_sel)); % select the indices of remaining vectors in the subset of ind_sel vectors 110 112 end 111 113 % if no value exceeds threshold, the result is recorded … … 116 118 y_dist=(Coord(ind_sel,2)CentreY(isub))/y_width;% relative ydistance to the retangle centre 117 119 weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge 120 %weight=1;% case for r1129 and before 118 121 U_smooth(ind_sel)=U_smooth(ind_sel)+weight.*U_smooth_sub; 119 122 V_smooth(ind_sel)=V_smooth(ind_sel)+weight.*V_smooth_sub; … … 138 141 y_dist=(Coord(ind_sel(ind_ind_sel),2)CentreY(isub))/y_width;% relative ydistance to the retangle centre 139 142 weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge 143 %weight=1; 140 144 U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+weight.*U_smooth_sub; 141 145 V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+weight.*V_smooth_sub; … … 166 170 U_smooth=U_smooth./nb_select;% take the average at the intersection of several subdomains 167 171 V_smooth=V_smooth./nb_select; 168 U_smooth(FF==4)=U(FF==4);% set to the initial values the eliminated vectors (flagged as false) 169 V_smooth(FF==4)=V(FF==4); 172 173 %eliminate the vectors with diff>threshold not yet eliminated 174 if exist('Threshold','var')&&~isempty(Threshold) 175 UDiff=U_smoothU;% difference between interpolated U component and initial value 176 VDiff=V_smoothV;% difference between interpolated V component and initial value 177 NormDiff=UDiff.*UDiff+VDiff.*VDiff;% Square of difference norm 178 FF(NormDiff>Threshold)=true;%put FF value to 1 to identify the criterium of elimmination 179 end 180 181 U_smooth(FF)=U(FF);% set to the initial values the eliminated vectors (flagged as false) 182 V_smooth(FF)=V(FF); 170 183 fill=zeros(NbCoord+1,NbCoord,size(SubRange,3)); %matrix of zeros to complement the matrix Data.Civ1_Coord_tps (conveninent for file storage) 171 184 Coord_tps=cat(1,Coord_tps,fill);
Note: See TracChangeset
for help on using the changeset viewer.