femwind/prolongation.m

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