femwind/midpoints2corners.m

   1 function C=midpoints2corners(M)
   2 % extrapolation formula
   3 extra = @(a,b) 2*a-b;
   4
   5 n=size(M);
   6 d=length(n);
   7 C=zeros(n+1);
   8 if d==1
   9     % inner-corners
  10     C(2:end-1)=avg1d(M);
  11     % corners
  12     C(1)=extra(C(2),C(3));
  13     C(end)=extra(C(end-1),C(end-2));
  14 elseif d==2
  15     % inner-corners
  16     C(2:end-1,2:end-1)=avg2d(M);
  17     % edge-corners
  18     C(1,2:end-1)=extra(C(2,2:end-1),C(3,2:end-1));
  19     C(end,2:end-1)=extra(C(end-1,2:end-1),C(end-2,2:end-1));
  20     C(2:end-1,1)=extra(C(2:end-1,2),C(2:end-1,3));
  21     C(2:end-1,end)=extra(C(2:end-1,end-1),C(2:end-1,end-2));
  22     % corners
  23     C(1,1)=extra(C(2,2),C(3,3));
  24     C(1,end)=extra(C(2,end-1),C(3,end-2));
  25     C(end,1)=extra(C(end-1,2),C(end-2,3));
  26     C(end,end)=extra(C(end-1,end-1),C(end-2,end-2));
  27 elseif d==3
  28     % inner-corners
  29     C(2:end-1,2:end-1,2:end-1)=avg3d(M);
  30     % face-corners
  31     C(1,2:end-1,2:end-1)=extra(C(2,2:end-1,2:end-1),C(3,2:end-1,2:end-1));
  32     C(end,2:end-1,2:end-1)=extra(C(end-1,2:end-1,2:end-1),C(end-2,2:end-1,2:end-1));
  33     C(2:end-1,1,2:end-1)=extra(C(2:end-1,2,2:end-1),C(2:end-1,3,2:end-1));
  34     C(2:end-1,end,2:end-1)=extra(C(2:end-1,end-1,2:end-1),C(2:end-1,end-2,2:end-1));
  35     C(2:end-1,2:end-1,1)=extra(C(2:end-1,2:end-1,2),C(2:end-1,2:end-1,3));
  36     C(2:end-1,2:end-1,end)=extra(C(2:end-1,2:end-1,end-1),C(2:end-1,2:end-1,end-2));
  37     % edge-corners
  38     C(1,1,2:end-1)=extra(C(2,2,2:end-1),C(3,3,2:end-1));
  39     C(end,1,2:end-1)=extra(C(end-1,2,2:end-1),C(end-2,3,2:end-1));
  40     C(1,end,2:end-1)=extra(C(2,end-1,2:end-1),C(3,end-2,2:end-1));
  41     C(end,end,2:end-1)=extra(C(end-1,end-1,2:end-1),C(end-2,end-2,2:end-1));
  42     C(1,2:end-1,1)=extra(C(2,2:end-1,2),C(3,2:end-1,3));
  43     C(end,2:end-1,1)=extra(C(end-1,2:end-1,2),C(end-2,2:end-1,3));
  44     C(1,2:end-1,end)=extra(C(2,2:end-1,end-1),C(3,2:end-1,end-2));
  45     C(end,2:end-1,end)=extra(C(end-1,2:end-1,end-1),C(end-2,2:end-1,end-2));
  46     C(2:end-1,1,1)=extra(C(2:end-1,2,2),C(2:end-1,3,3));
  47     C(2:end-1,end,1)=extra(C(2:end-1,end-1,2),C(2:end-1,end-2,3));
  48     C(2:end-1,1,end)=extra(C(2:end-1,2,end-1),C(2:end-1,3,end-2));
  49     C(2:end-1,end,end)=extra(C(2:end-1,end-1,end-1),C(2:end-1,end-2,end-2));
  50     % corners
  51     C(1,1,1)=extra(C(2,2,2),C(3,3,3));
  52     C(end,1,1)=extra(C(end-1,2,2),C(end-2,3,3));
  53     C(1,end,1)=extra(C(2,end-1,2),C(3,end-2,3));
  54     C(1,1,end)=extra(C(2,2,end-1),C(3,3,end-2));
  55     C(1,end,end)=extra(C(2,end-1,end-1),C(3,end-2,end-2));
  56     C(end,1,end)=extra(C(end-1,2,end-1),C(end-2,3,end-2));
  57     C(end,end,1)=extra(C(end-1,end-1,2),C(end-2,end-2,3));
  58     C(end,end,end)=extra(C(end-1,end-1,end-1),C(end-2,end-2,end-2));
  59 else
  60     error(['Not implemented number of dimensions ',d]);
  61 end
  62
  63
  64 function A=avg1d(D)
  65     A=( D(1:end-1) + D(2:end) ) / 2;
  66 end
  67
  68 function A=avg2d(D)
  69     A=( D(1:end-1,1:end-1) + D(1:end-1,2:end) + ...
  70         D(2:end,1:end-1) + D(2:end,2:end) ) / 4;
  71 end
  72
  73 function A=avg3d(D)
  74     A=( D(1:end-1,1:end-1,1:end-1) + D(1:end-1,1:end-1,2:end) + ...
  75         D(1:end-1,2:end,1:end-1) + D(1:end-1,2:end,2:end) + ...
  76         D(2:end,1:end-1,1:end-1) + D(2:end,1:end-1,1:end-1) + ...
  77         D(1:end-1,1:end-1,1:end-1) + D(1:end-1,1:end-1,1:end-1) ) / 8;
  78 end
  79 end