interpolate.m

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

%========================================================================
%========================================================================
%INTERPOLATE Transfer correction from the coarse grid to the current grid.
%========================================================================
%      This is a beauty datas of interpolation in MG made by zhiyong 
%  It provides several interpolation methods which can be used as transfer
%  operator of error in MG algorithm.
%========================================================================
%                Zhiyong 2008 starting this work 
%            through considering kinds of interpolation
% our mesh ordered by from south to north ,from left to right
%                 ---③--⑥--⑨--
%                    |   |   |
%                 ---②--⑤--⑧--
%                    |   |   |
%                  --①--④--⑦--
function correct = interpolate(level, u_c_out,rr)

include_flags
include_globals 
extract_globals
IX_c = 1:(nx_c+1); IY_c = 1:(ny_c+1);
IX_f = 1:(nx_f+1); IY_f = 1:(ny_f+1);


if interp_flag == LINEAR

   COARSE(IX_c,IY_c) = reshape(u_c_out,nx_c+1,ny_c+1);
   FINE = interp2(X_c,Y_c,COARSE,X_f,Y_f);
   correct = reshape(FINE(IX_f,IY_f), N_f, 1);
%===================两种差分方法结合解局部方程======================
elseif interp_flag == FULL_LOCAL_EQUATION
      nx_c=nx_c+1;
      ny_c=nx_c;
      nx_f=nx_f+1;
      ny_f=nx_f;
      h=1 / (nx_f-1);  % mesh size on current level
     right = reshape(rr,nx_f,ny_f);
     coarse(1:nx_c,1:ny_c) = reshape(u_c_out,nx_c,ny_c);
     fine(1:nx_f,1:ny_f)=zeros(nx_f,ny_f);
     for j=1:nx_c
         for k=1:ny_c
             fine(2*k-1,2*j-1)=coarse(k,j);
         end
     end        %============ points in coarse mesh ==============
     for j=1:nx_c-1
         for k=1:ny_c-1
             fine(2*k,2*j)=0.25*(coarse(k,j)+coarse(k+1,j)+...
                 coarse(k,j+1)+coarse(k+1,j+1)+h*h*right(2*k,2*j));
         end
     end  %========== points isnot in coarse mesh==================
     
     for j=1:nx_c-2
         for k=1:ny_c-1
             fine(2*k,2*j+1)=0.25*(coarse(k+1,j+1)+coarse(k,j+1)+...
                 fine(2*k,2*j)+fine(2*k,2*j+2)+h*h*right(2*k,2*j+1));
         end
     end%  ???points between two  northsouth coarse points=========
    for j=1:nx_c-1
         for k=1:ny_c-2
             fine(2*k+1,2*j)=0.25*(coarse(k+1,j+1)+coarse(k+1,j)+...
                 fine(2*k,2*j)+fine(2*k+2,2*j)+h*h*right(2*k+1,2*j));
         end
     end% ???points between two we direction coarse points=========
     
     for j=1
         for k=1:ny_c-1
             fine(2*k,2*j-1)=0.5*(coarse(k,j)+coarse(k+1,j));
         end
     end  % 西===========================================
     
     for j=nx_c
         for k=1:ny_c-1
             fine(2*k,2*j-1)=0.5*(coarse(k,j)+coarse(k+1,j));
         end
     end   % 东=====================================
     for k=1
         for j=1:nx_c-1
             fine(2*k-1,2*j)=0.5*(coarse(k,j)+coarse(k,j+1));
         end
     end  %北=====================================
   for k=ny_c
         for j=1:nx_c-1
             fine(2*k-1,2*j)=0.5*(coarse(k,j)+coarse(k,j+1));
         end
     end   % 南===================================
     
     correct = reshape(fine,N_f,1);
end