femwind/check_nonzeros.m

   1 function check_nonzeros(level,Kc,Xc,P,K,X)
   2 % check_nonzeros(Kc,Xc,P,K,X)
   3 % checking nonzeros structure consistent with hexa grid
   4     disp('check_nonzeros: check if structure consistent with hexa grid')
   5     nc=size(Xc{1});
   6     nnc=prod(nc);
   7     fprintf('level %g size %g %g %g stiffness matrix size %g nonzeros %g density %g%%\n',...
   8         level,nc,length(Kc),nnz(Kc),100*nnz(Kc)/prod(size(Kc)))
   9     if exist('P','var')
  10         fprintf('prolongation matrix size %g %g nonzeros %g density %g%%\n',...
  11             size(P),nnz(P),100*nnz(P)/prod(size(P)))
  12         fprintf('prolongation matrix max row nonzeros %g max column nonzeros %g\n',...
  13             max(full(sum(P~=0,2))),max(full(sum(P~=0,1))))
  14     end
  15     if any(size(Kc)~=nnc)
  16         error('Kc size inconsistent with mesh Xc')
  17     end
  18     if any(any(sign(Kc)-sign(Kc)'))
  19         error('K does not have symmetric structure')
  20     end
  21     [ii,jj,aij]=find(Kc);
  22     [i1,i2,i3]=ind2sub(nc,ii);
  23     [j1,j2,j3]=ind2sub(nc,jj);
  24     d=max([abs(i1-j1),abs(i2-j2),abs(i3-j3)],[],2);
  25     f=find(d>1);
  26     if f
  27         for xj=f(:)'
  28             i = ii(xj);
  29             j = jj(xj);
  30             fprintf('%i mesh node %i %i %i and %i mesh node %i %i %i distance %i value %g\n',...
  31                 i,i1(xj),i2(xj),i3(xj),...
  32                 j,j1(xj),j2(xj),j3(xj),d(xj),aij(xj))
  33             if exist('P','var') & exist('X','var')
  34                 n = size(X{1});
  35                 fi = find(P(:,i))';
  36                 [fi1,fi2,fi3]=ind2sub(n,fi);
  37                 [i1,i2,i3]=ind2sub(n,i);
  38                 fj = find(P(:,j))';
  39                 [fj1,fj2,fj3]=ind2sub(n,fj);
  40                 [j1,j2,j3]=ind2sub(n,j);
  41                 fprintf('coarse %i at %i %i %i interpolates to\n',i,i1,i2,i3)
  42                 for k=1:length(fi)
  43                     fprintf('fine   %i at %i %i %i\n',fi(k),fi1(k),fi2(k),fi3(k))
  44                 end
  45                 fprintf('coarse %i at %i %i %i interpolates to\n',j,j1,j2,j3)
  46                 for k=1:length(fj)
  47                     fprintf('fine   %i at %i %i %i\n',fj(k),fj1(k),fj2(k),fj3(k))
  48                 end
  49                 disp('connected by nonzeros of fine matrix')
  50                 [fii,fjj,fval]=find(K(fi,fj));
  51                 sub=sparse(fi(fii),fj(fjj),fval);
  52                 disp(sub)
  53                 error('K inconsistent with hexa structure nonzeros at distance more than 1')
  54             end
  55         end
  56     end
  57 end
  58