interpolate.asv

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

%========================================================================
%========================================================================
%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);
h=1 / (nx_f+1);  % mesh size on current level
[t1,t2] = meshgrid([h:h:(1-h)]);
%UTRUE = sin(2*pi*t1) .* sin(3*pi*t2);
%right= (h*h) * 13*pi*pi*UTRUE;   % right of equation 
%w2=1;w1=1;
%UTRUE = sin(2*pi*w1*t1) .* sin(3*pi*w2*t2);   % solution of equation 
%right= (h*h) * (4*w1*w1+9*w2*w2)*pi*pi*UTRUE;   % right of equation


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 == CUBIC

   %COARSE(IX_c,IY_c) = reshape(u_c_out,nx_c,ny_c);
   %FINE = interp2(X_c,Y_c,COARSE,X_f,Y_f,'cubic');
   %correct = reshape(FINE(IX_f,IY_f), N_f, 1);

elseif interp_flag == EXPLICIT_BILINEAR

   eval(['PROLONG = ARRAY',num2str(level), ';']);
   correct = PROLONG * u_c_out;
%***********adding my new interpolation in MG************************
%===================方法一：两种差分方法结合解局部方程======================
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=ny_f
         for j=2:nx_f-1
             fine(k,j)=0.5*fine(k-1,j);
         end
     end   %北=====================================
     for k=1
         for j=2:nx_f-1
             fine(k,j)=0.5*fine(k+1,j);
         end
     end  % 南===================================
     
     correct = reshape(fine,N_f,1);
end