restrict.m

五点差分型多重网格方法：各种插值算子的比较）

%RESTRICT Transfer residual from the current grid to the next coarser grid.
%
%       RHS_C = RESTRICT(LEVEL,R) uses the restriction scheme defined by 
%       "restrict_flag" to transfer the vector R on the current level LEVEL 
%       to the vector RHS_C on the next coarser level LEVEL+1.
% =======================================================================
%  Here I will provide several methods according to our interpolations.
%  our mesh ordered by from south to north ,from left to right
%                 ---③--⑥--⑨--
%                    |   |   |
%                 ---②--⑤--⑧--
%                    |   |   |
%                  --①--④--⑦--
%①、③、⑦、⑨ are coarse points, we need to obtain their values.
% Zhiyong will finish it before 20008.10.28
function rhs_c = restrict(level,r)

include_globals 
extract_globals
include_flags
%--------------others'design in starting---------------------
    nx_f=nx_f+1;
    ny_f=nx_f;
    nx_c=nx_c+1;
    ny_c=nx_c;
%--------------------------------------------------------------    
   nx0_f = nx_f+3;
   ny0_f = ny_f+3;
   N0_f = nx0_f*ny0_f;
   dx=1;
   dy=nx0_f;
   r0 = zeros(N0_f,1);
   for iy=1:ny_f
   for ix=1:nx_f
       r0(nx0_f+1 + ix + nx0_f*(iy-1)) = r(ix+nx_f*(iy-1));
   end
   end
   I = zeros(N_c,1);
   for iy=1:ny_c
   for ix=1:nx_c
       I(ix + nx_c*(iy-1)) = 2*ix + 2*iy*nx0_f + 1;
   end
   end
if restrict_flag == INJECTION

   rhs_c = r0(I);
   rhs_c = 4*rhs_c;
elseif restrict_flag == HALF_WEIGHTING

   rhs_c = .5*r0(I) + ...
          .125*(r0(I+dx) + r0(I-dx) + r0(I+dy) + r0(I-dy));
      
   rhs_c = 4*rhs_c;
elseif restrict_flag == FULL_WEIGHTING

   rhs_c = .25*r0(I) + ...
          .125*(r0(I+dx) + r0(I-dx) + r0(I+dy) + r0(I-dy)) + ...
          .0625*(r0(I+dx+dy) + r0(I-dx+dy) + r0(I+dx-dy) + r0(I-dx-dy));
   rhs_c = 4*rhs_c;   
%=================方法一：直接作一个映射从细网格到粗网格=====================
elseif restrict_flag == ZHIJIE
       
     fine(1:nx_f,1:ny_f)=reshape(r,nx_f,ny_f); 
     coarse(1:nx_c,1:ny_c)=zeros(nx_c,ny_c);
     for j=1:nx_c-1
        for k=1:ny_c-1
           coarse(k,j)=fine(2*k,2*j);
        end
     end
     rhs_c=reshape(coarse,(nx_c)*(ny_c),1);
     rhs_c = 4*rhs_c;
end