femwind/sparse_assembly.m

   1 function [K,F,W]=sparse_assembly(A,X,u0,lambda,params)
   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; if empty, do not compute K and F
   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 % KK = sparse([],[],[],nn,nn,nz);
  24 F = zeros(nn,1);
  25 W = {zeros(n-1),zeros(n-1),zeros(n-1)};
  26 nzelem=prod(n-1)*64;     % total nonzeros before assembly
  27 ii=zeros(nzelem,1);
  28 jj=zeros(nzelem,1);
  29 aa=zeros(nzelem,1);
  30 kk=0;
  31
  32 % create matrices
  33 tstart=tic;
  34 Xloc = zeros(3,8);
  35 kglo=zeros(1,8);
  36 m = n-1;
  37 done_last=0;
  38 disp('accumulating global stiffness matrix entries')
  39
  40 for i3=1:m(3)
  41     for i2=1:m(2)
  42         for i1=1:m(1)  % loop over elements
  43             % now in element (i1,i2,i3)
  44             for j3=0:1 % loop over corners of the element
  45                 for j2=0:1
  46                     for j1=0:1
  47                         kloc=1+j1+2*(j2+2*j3);  % local index of the node in element (i1, i2. i3)
  48                         % kloc1=sub2ind([2,2,2],j1+1,j2+1,j3+1);  if kloc1~=kloc, error('kloc'),end
  49                         k1 = i1+j1; k2 = i2+j2; k3 = i3+j3; %  position of the node in the global grid
  50                         kglo(kloc)=k1+n(1)*((k2-1)+n(2)*(k3-1)); % global index
  51                         % kglo1=sub2ind(n,i1+j1,i2+j2,i3+j3); if kglo1 ~= kglo(kloc), error('kglo'), end
  52                         for i=1:3
  53                             Xloc(i,kloc)=X{i}(i1+j1,i2+j2,i3+j3); % node coordinate i
  54                         end
  55                     end
  56                 end
  57             end
  58             if ~isempty(u0)
  59                 u0loc=[u0{1}(i1,i2,i3),u0{2}(i1,i2,i3),u0{3}(i1,i2,i3)]';
  60             else
  61                 u0loc=[];
  62             end
  63             [Kloc,Floc,Jg]=hexa(A,Xloc,u0loc); % create local matrix and rhs
  64             F(kglo)=F(kglo)+Floc; % assemble to global rhs
  65             % KK(kglo,kglo)=KK(kglo,kglo)+Kloc; % assemble to global matrix
  66             % instead, accumulate contributions to the global matrix
  67             [ix,jx]=ndgrid(kglo,kglo);
  68             nzloc=prod(size(Kloc));
  69             ii(kk+1:kk+nzloc)=ix(:);
  70             jj(kk+1:kk+nzloc)=jx(:);
  71             aa(kk+1:kk+nzloc)=Kloc(:);
  72             kk=kk+nzloc;
  73             if ~isempty(lambda)
  74                 grad = lambda(kglo)'*Jg;  % grad lambda
  75                 grad = grad/A;
  76                 for i=1:3
  77                     W{i}(i1,i2,i3)=grad(i);
  78                 end
  79             end
  80         end
  81         % done = 100*((i3-1)*m(2)+i2)/(m(3)*m(2));
  82         % done = round(done);
  83         % if done>done_last+5, fprintf(' %g%% ',done), done_last=done; end
  84     end
  85 end
  86
  87 if ~isempty(u0)
  88     for i=1:3
  89         W{i}=u0{i}+W{i};
  90     end
  91 end
  92 % fprintf('\n')
  93 if kk~=nzelem, error('wrong element nonzeros'),end
  94 disp('building sparse global stiffness matrix')
  95 K = sparse(ii,jj,aa,nn,nn); % err_K=big(K-KK)
  96 nn=prod(n);
  97 % checks
  98 if length(F)>nn, size(F),error('size has grown by an index out of bounds'),end
  99 if exist('params','var') && isfield(params,'lambda')
 100     check_nonzeros(params.levels,K,X)
 101 end
 102 check_symmetry(K,'K',eps)
 103 K = (K+K')/2; % symmetrize
 104 end