femwind/coarsening_P.m

   1 function P=coarsening_P(hzc,icl3,X,params)
   2 % [P,X_coarse]=coarsening_P(icl,X)
   3 % Build the prolongation matrix
   4 % In:
   5 %   hzc     coarsening factor in the horizontal directions 1 and 2
   6 %   icl     indices of coarse grid in the 3 directions
   7 %   X       grid coordinates
   8 %   params
   9
  10     n = size(X{1});   % grid size
  11     nn = prod(n);     % number of grid points
  12     [icl1,icl2]=hzc2icl(hzc,n);
  13     nc = [length(icl1),length(icl2),length(icl3)];
  14     nnc = prod(nc);  % number of coarse grid points
  15
  16     disp('building the prolongation')
  17     nncz=nnc*27; % estimate number of nonzeros in coarse matrix
  18     ia=zeros(nncz,1); % preallocate arrays for P matrix entries
  19     ja=ia;aa=ia;
  20     % P = spalloc(nn,nnc,nnc*27);
  21     k=0;
  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 ic1=1:nc(1)        % horizontal loops over coarse points
  27             if1=icl1(ic1);  % fine grid index of the coarse point
  28             [ifs1,ife1]=ifse(ic1,icl1,nc(1)); % get start and end of support
  29             % fprintf('coarse x1 %g at %g contributes to %g : %g\n',ic1,if1,ifs1,ife1)
  30             for ic2=1:nc(2)
  31                 if2=icl2(ic2);  % fine grid index of the coarse point
  32                 [ifs2,ife2]=ifse(ic2,icl2,nc(2)); % get start and end of support
  33                 % fprintf('coarse x2 %g at %g contributes to %g : %g\n',ic2,if2,ifs2,ife2)
  34                 ixc = sub2ind(nc,ic1,ic2,ic3); % index into coarse matrix
  35                 % loop over fine points coarse point ic1 ic2 ic3
  36                 % 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                     for i2=ifs2:ife2
  42                         for i3=ifs3:ife3
  43                             ix=sub2ind(n,i1,i2,i3);
  44                             if params.P_by_x
  45                                 % should adapt from X
  46                                 if i1>if1
  47                                     q1=(X{1}(i1,i2,i3)-X{1}(ife1+1,i2,i3))/(X{1}(if1,i2,i3)-X{1}(ife1+1,i2,i3));
  48                                 elseif i1<if1
  49                                     q1=(X{1}(i1,i2,i3)-X{1}(ifs1-1,i2,i3))/(X{1}(if1,i2,i3)-X{1}(ifs1-1,i2,i3));
  50                                 else
  51                                     q1=1;
  52                                 end
  53                                 if i2>if2
  54                                     q2=(X{2}(i1,i2,i3)-X{2}(i1,ife2+1,i3))/(X{2}(i1,if2,i3)-X{2}(i1,ife2+1,i3));
  55                                 elseif i2<if2
  56                                     q2=(X{2}(i1,i2,i3)-X{2}(i1,ifs2-1,i3))/(X{2}(i1,if2,i3)-X{2}(i1,ifs2-1,i3));
  57                                 else
  58                                     q2=1;
  59                                 end
  60                                 if i3>if3
  61                                     q3=(X{3}(i1,i2,i3)-X{3}(i1,i2,ife3+1))/(X{3}(i1,i2,if3)-X{3}(i1,i2,ife3+1));
  62                                 elseif i3<if3
  63                                     q3=(X{3}(i1,i2,i3)-X{3}(i1,i2,ifs3-1))/(X{3}(i1,i2,if3)-X{3}(i1,i2,ifs3-1));
  64                                 else
  65                                     q3=1;
  66                                 end
  67                                 val=q1*q2*q3;
  68                             else
  69                                 val=(1-0.5*abs(i1-if1))...
  70                                     *(1-0.5*abs(i2-if2))...
  71                                     *(1-0.5*abs(i3-if3)); % horizontal
  72                             end
  73                             k=k+1;
  74                             ia(k)=ix;
  75                             ja(k)=ixc;
  76                             aa(k)=val;
  77                         end
  78                     end
  79                 end
  80             end
  81         end
  82     end
  83     P = sparse(ia(1:k),ja(1:k),aa(1:k),nn,nnc);
  84 end
  85
  86
  87 function [ifs,ife]=ifse(ic,icl,ncn)
  88     if ic>1
  89         ifs=icl(ic-1)+1; % from above previous coarse
  90     else
  91         ifs=icl(ic);     % itself, there is no previous fine layer
  92     end
  93     if ic<ncn
  94         ife=icl(ic+1)-1; % up to under next layer
  95     else
  96         ife=icl(ic);     % itself, there is no next layer
  97     end
  98 end