femwind/sparse_assembly.m

   1 function [K,F,W]=sparse_assembly(A,X,u0,lambda)
   2 % in:
   3 %  A  penalty coefficients matrix, size 3x3, s.p.d.
   4 %  X  cell array with 3d arrays x y z coordinates of mesh vertices
   5 %  u0 cell array with 3d arrays initial wind vector at centers in x y z directions
   6 %  lambda scalar field on node
   7 % out:
   8 %  K  stiffness matrix, sparse
   9 %  F  load vector
  10 %  W  gradient of lambda
  11
  12 % do to: now computing all three K,F,W, add a switch to compute only one
  13 % and save time
  14
  15 fprintf('sparse_assembly:')
  16
  17 d = size(X,2);    % dimensions
  18 n = size(X{1});   % mesh size
  19 nn = prod(n);     % total nodes
  20 nz = nn*27;        % estimated nonzeros
  21
  22 % initialize matrices
  23 K = sparse([],[],[],nn,nn,nz);
  24 F = zeros(nn,1);
  25 W = {zeros(n-1),zeros(n-1),zeros(n-1)};
  26
  27 % create matrices
  28 tstart=tic;
  29 Xloc = zeros(3,8);
  30 kglo=zeros(1,8);
  31 m = n-1;
  32 done_last=0;
  33 for i3=1:m(3)
  34     for i2=1:m(2)
  35         for i1=1:m(1)  % loop over elements
  36             % now in element (i1,i2,i3)
  37             for j3=0:1 % loop over corners of the element
  38                 for j2=0:1
  39                     for j1=0:1
  40                         kloc=sub2ind([2,2,2],j1+1,j2+1,j3+1); % local index
  41                         % kloc=1+j1+2*(j2+2*j3);
  42                         kglo(kloc)=sub2ind(n,i1+j1,i2+j2,i3+j3); % global index
  43                         for i=1:3
  44                             Xloc(i,kloc)=X{i}(i1+j1,i2+j2,i3+j3); % node coordinate i
  45                         end
  46                     end
  47                 end
  48             end
  49             u0loc=[u0{1}(i1,i2,i3),u0{2}(i1,i2,i3),u0{3}(i1,i2,i3)]';
  50             [Kloc,Floc,Jg]=hexa(A,Xloc,u0loc); % create local matrix and rhs
  51             K(kglo,kglo)=K(kglo,kglo)+Kloc; % assemble to global matrix
  52             F(kglo)=F(kglo)+Floc; % assemble to global rhs
  53             grad = lambda(kglo)'*Jg;  % grad lambda
  54             grad = grad/A;
  55             for i=1:3
  56                 W{i}(i1,i2,i3)=grad(i);
  57             end
  58         end
  59         done = 100*((i3-1)*m(2)+i2)/(m(3)*m(2));
  60         done = round(done);
  61         if done>done_last+5
  62             fprintf(' %g%% ',done)
  63             done_last=done;
  64         end
  65     end
  66 end
  67 for i=1:3
  68     W{i}=u0{i}+W{i};
  69 end
  70 fprintf('\n')
  71 nn=prod(n);
  72 nz=nnz(K);
  73 fprintf('stiffness matrix size %g nonzeros %g density %g%%\n',nn,nz,100*nz/nn^2)
  74 % checks
  75 if length(F)>nn, size(F),error('size has grown by an index out of bounds'),end
  76 check_symmetry(K,'K',eps)
  77 K = (K+K')/2; % symmetrize
  78 end