quicwind/restriction_3d.m

   1 function y=restriction_3d(x)
   2     % in:
   3     %    x   array
   4     % out
   5     %    y   bilinear average on twice coarser grid
   6     % this is transpose of prolongation with weighting
   7
   8     % to test:
   9     % x=magic(2); restriction_3d(x)
  10
  11     % map vector to 3d grid with zero boundary
  12     % average from neighbors with the same weights as prolongation
  13     % scaled to sum one
  14     if any(mod(size(x),2))==0
  15         error('restriction_3d: input dimensions must be odd')
  16     end
  17     tw=1/(1+6*1/2+12*1/4+8*1/8);
  18     i1=0; i2=0; i3=0;
  19     y = tw*x(2+i1:2:end-1+i1,2+i2:2:end-1+i2,2+i3:2:end-1+i3);
  20     for i1=-1:1
  21         for i2=-1:1
  22             for i3=-1:1
  23                 if (i1 ~=0) || (i2 ~=0) || (i3 ~=0)
  24                     w = tw/((1+abs(i1))*(1+abs(i2))*(1+abs(i3)));
  25                     y = y + w*x(2+i1:2:end-1+i1,2+i2:2:end-1+i2,2+i3:2:end-1+i3);
  26                 end
  27             end
  28         end
  29     end
  30 end