femwind/restriction.m

   1 function uc=restriction(u,hzc,icl3,X,params)
   2 % [P,X_coarse]=restriction(icl,X)
   3 % Multiply by prolongation matrix transpose
   4 % In:
   5 %   u       fine grid vector
   6 %   icl     cell array size 3, indices of coarse grid in the 3 directions
   7 %   X       grid coordinates
   8 %   params
   9 %
  10 % Out:
  11 %   uc      coarse grid vector
  12
  13     n = size(X{1});   % grid size
  14     nn = prod(n);     % number of grid points
  15     [icl1,icl2]=hzc2icl(hzc,n);
  16     nc = [length(icl1),length(icl2),length(icl3)];
  17     nnc = prod(nc);  % number of coarse grid points
  18
  19     uc = zeros(nc);  % allocate
  20     disp('restriction')
  21     for ic3=1:nc(3)           % loop over coarse layers
  22         if3=icl3(ic3);         % the fine number of the coarse layer
  23         [ifs3,ife3]=ifse(ic3,icl3,nc(3)); % get start and end of support
  24         % fprintf('vertical coarse layer %g at %g contributes to layers %g : %g\n',ic3,if3,ifs3,ife3)
  25         for ic1=1:nc(1)        % horizontal loops over coarse points
  26             if1=icl1(ic1);  % fine grid index of the coarse point
  27             [ifs1,ife1]=ifse(ic1,icl1,nc(1)); % get start and end of support
  28             % fprintf('coarse x1 %g at %g contributes to %g : %g\n',ic1,if1,ifs1,ife1)
  29             for ic2=1:nc(2)
  30                 if2=icl2(ic2);  % fine grid index of the coarse point
  31                 [ifs2,ife2]=ifse(ic2,icl2,nc(2)); % get start and end of support
  32                 % fprintf('coarse x2 %g at %g contributes to %g : %g\n',ic2,if2,ifs2,ife2)
  33                 % loop over fine points coarse point ic1 ic2 ic3 contributes to
  34                 % coarse point ic1 ic2 ic3 is if1 if2 if3 on the fine grid
  35                 % interpolating from between to (i1 i2 i3) from if1 if2 if3
  36                 % and the next coarse
  37                 for i1=ifs1:ife1
  38                     for i2=ifs2:ife2
  39                         for i3=ifs3:ife3
  40                                 if i1>if1
  41                                     q1=(X{1}(i1,i2,i3)-X{1}(ife1+1,i2,i3))/(X{1}(if1,i2,i3)-X{1}(ife1+1,i2,i3));
  42                                 elseif i1<if1
  43                                     q1=(X{1}(i1,i2,i3)-X{1}(ifs1-1,i2,i3))/(X{1}(if1,i2,i3)-X{1}(ifs1-1,i2,i3));
  44                                 else
  45                                     q1=1;
  46                                 end
  47                                 if i2>if2
  48                                     q2=(X{2}(i1,i2,i3)-X{2}(i1,ife2+1,i3))/(X{2}(i1,if2,i3)-X{2}(i1,ife2+1,i3));
  49                                 elseif i2<if2
  50                                     q2=(X{2}(i1,i2,i3)-X{2}(i1,ifs2-1,i3))/(X{2}(i1,if2,i3)-X{2}(i1,ifs2-1,i3));
  51                                 else
  52                                     q2=1;
  53                                 end
  54                                 if i3>if3
  55                                     q3=(X{3}(i1,i2,i3)-X{3}(i1,i2,ife3+1))/(X{3}(i1,i2,if3)-X{3}(i1,i2,ife3+1));
  56                                 elseif i3<if3
  57                                     q3=(X{3}(i1,i2,i3)-X{3}(i1,i2,ifs3-1))/(X{3}(i1,i2,if3)-X{3}(i1,i2,ifs3-1));
  58                                 else
  59                                     q3=1;
  60                                 end
  61                                 uc(ic1,ic2,ic3) = uc(ic1,ic2,ic3) + q1*q2*q3*u(i1,i2,i3);
  62                         end
  63                     end
  64                 end
  65             end
  66         end
  67     end
  68 end
  69
  70
  71 function [ifs,ife]=ifse(ic,icl,ncn)
  72     if ic>1
  73         ifs=icl(ic-1)+1; % from above previous coarse
  74     else
  75         ifs=icl(ic);     % itself, there is no previous fine layer
  76     end
  77     if ic<ncn
  78         ife=icl(ic+1)-1; % up to under next layer
  79     else
  80         ife=icl(ic);     % itself, there is no next layer
  81     end
  82 end