quicwind/prolongation_3d.m

   1 function y=prolongation_3d(x)
   2     % in:
   3     %    x   3d array
   4     %    nn  shape of x
   5     % out
   6     %    y   bilinear interpolation wrapped by zeros
   7
   8     % to test:
   9     % x=magic(2); prolongation_3d(x)
  10
  11     % map vector to 3d grid with zero boundary
  12     nn=size(x);
  13     xx=zeros(nn+2);
  14     xx(2:nn(1)+1,2:nn(2)+1,2:nn(3)+1)=reshape(x,nn);
  15     % allocate output
  16     y=zeros(2*nn+1);
  17     % copy values on coarse points
  18     y(2:2:end-1,2:2:end-1,2:2:end-1)=xx(2:end-1,2:end-1,2:end-1);
  19
  20     % averages to coarse edge midpoints in direction 1
  21     y(1:2:end,2:2:end-1,2:2:end-1)=0.5*(xx(1:end-1,2:end-1,2:end-1)+xx(2:end,2:end-1,2:end-1));
  22     % averages to coarse edge midpoints in direction 2
  23     y(2:2:end-1,1:2:end,2:2:end-1)=0.5*(xx(2:end-1,1:end-1,2:end-1)+xx(2:end-1,2:end,2:end-1));
  24     % averages to coarse edge midpoints in direction 3
  25     y(2:2:end-1,2:2:end-1,1:2:end)=0.5*(xx(2:end-1,2:end-1,1:end-1)+xx(2:end-1,2:end-1,2:end));
  26
  27     % averages to coarse cell face centers in directions 1 and 2
  28     y(1:2:end,1:2:end,2:2:end-1)=0.25*(xx(1:end-1,1:end-1,2:end-1)+xx(2:end,1:end-1,2:end-1)+...
  29                              xx(1:end-1,2:end,2:end-1) + xx(2:end,2:end,2:end-1));
  30     % averages to coarse cell face centers in directions 2 and 3
  31     y(1:2:end,2:2:end-1,1:2:end)=0.25*(xx(1:end-1,2:end-1,1:end-1)+xx(2:end,2:end-1,1:end-1)+...
  32                              xx(1:end-1,2:end-1,2:end) + xx(2:end,2:end-1,2:end));
  33     % averages to coarse cell face centers in directions 2 and 3
  34     y(2:2:end-1,1:2:end,1:2:end)=0.25*(xx(2:end-1,1:end-1,1:end-1)+xx(2:end-1,2:end,1:end-1)+...
  35                              xx(2:end-1,1:end-1,2:end) + xx(2:end-1,2:end,2:end));
  36
  37     % averages to coarse cell centers
  38     y(1:2:end,1:2:end,1:2:end)=0.125*(xx(1:end-1,1:end-1,1:end-1)+xx(2:end,1:end-1,1:end-1)+...
  39                              xx(1:end-1,2:end,1:end-1) + xx(1:end-1,1:end-1,2:end)+...
  40                              xx(2:end,2:end,1:end-1) + xx(2:end,1:end-1,2:end)+...
  41                              xx(1:end-1,2:end,2:end) + xx(2:end,2:end,2:end));
  42
  43 end