femwind/smoothing.m

   1 function x=smoothing(K,F,X,x,params)
   2 % x=smoothing(K,F,X,x)
   3 % one smoothing iteration
   4     n = size(X{1});
   5     it_type=sprintf('level %g smoothing',params.levels);
   6     switch params.smoothing
   7         case {'horizontal planes'}
   8             % disp('red-black vertical, planes horizontal')
   9             for rb3=1:2
  10                 for i3=rb3:2:n(3)
  11                         [planex,planey]=ndgrid(1:n(1),1:n(2));
  12                         planez = i3*ones(n(1),n(2));
  13                         % solving horizontal layer
  14                         ix = sub2ind(n,planex,planey,planez);
  15                         x(ix) = x(ix) - K(ix,ix)\(K(:,ix)'*x - F(ix));
  16                 end
  17             end
  18         case {'vertical lines'}
  19             % disp('red-black relaxation horizontal, lines vertical')
  20             for rb1=1:2
  21                 for rb2=1:2
  22                     for i1=rb1:2:n(1)
  23                         for i2=rb2:2:n(2)
  24                             % solving horizontal location i1 i2 and vertical line
  25                             ix = sub2ind(n,i1*onez,i2*onez,colz);
  26                             x(ix) = x(ix) - K(ix,ix)\(K(:,ix)'*x - F(ix));
  27                         end
  28                     end
  29                 end
  30             end
  31         case {'vertical sweeps'}
  32             % disp('red-black relaxation horizontal, down to up sweep vertical')
  33             x = vertical_sweeps(K,F,X,x);
  34         case '3D red-black'
  35             for rb1=1:2
  36                 for rb2=1:2
  37                     for rb3=1:2
  38                         for i1=rb1:2:n(1)
  39                             for i2=rb2:2:n(2)
  40                                 for i3=rb3:2:n(3)
  41                                     ix = sub2ind(n,i1,i2,i3);
  42                                     x(ix) = x(ix) - K(ix,ix)\(K(:,ix)'*x - F(ix));
  43                                 end
  44                             end
  45                         end
  46                     end
  47                 end
  48             end
  49         otherwise
  50             error(['smoothing ',params.smoothing,' unknown'])
  51     end
  52 end