femwind/w_assembly.m

   1 function W=w_assembly(A,X,u0,lambda,params)
   2
   3 % in:
   4 %  A  penalty coefficients matrix, size 3x3, s.p.d.
   5 %  X  cell array with 3d arrays x y z coordinates of mesh vertices
   6 %  u0 cell array with 3d arrays initial wind vector at centers in x y z directions
   7 %  lambda scalar field on node; if empty, do not compute K and F
   8 % out:
   9 %  K  stiffness matrix, sparse
  10 %  F  load vector
  11 %  W  gradient of lambda
  12
  13 % do to: now computing all three K,F,W, add a switch to compute only one
  14 % and save time
  15
  16 fprintf('w_assembly:')
  17
  18 n = size(X{1});   % mesh size
  19
  20 % initialize matrices
  21 W = {zeros(n-1),zeros(n-1),zeros(n-1)};
  22
  23 % create matrices
  24 Xloc = zeros(3,8);
  25 kglo=zeros(1,8);
  26 m = n-1;
  27 done_last=0;
  28 %disp('Constructing W')
  29 %disp('lambda array is')
  30 %disp(lambda)
  31 %disp('u0 from matlab is')
  32 %disp(u0{1})
  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=1+j1+2*(j2+2*j3);  % local index of the node in element (i1, i2. i3)
  41                         % kloc1=sub2ind([2,2,2],j1+1,j2+1,j3+1);  if kloc1~=kloc, error('kloc'),end
  42                         k1 = i1+j1; k2 = i2+j2; k3 = i3+j3; %  position of the node in the global grid
  43                         kglo(kloc)=k1+n(1)*((k2-1)+n(2)*(k3-1)); % global index
  44                         % kglo1=sub2ind(n,i1+j1,i2+j2,i3+j3); if kglo1 ~= kglo(kloc), error('kglo'), end
  45                         for i=1:3
  46                             Xloc(i,kloc)=X{i}(i1+j1,i2+j2,i3+j3); % node coordinate i
  47                         end
  48                     end
  49                 end
  50             end
  51             %disp('Xloc is')
  52             %disp(Xloc)
  53             if ~isempty(u0)
  54                 u0loc=[u0{1}(i1,i2,i3),u0{2}(i1,i2,i3),u0{3}(i1,i2,i3)]';
  55             else
  56                 u0loc=[];
  57             end
  58            %disp('u0loc is')
  59           % disp(u0loc(1))
  60             [~,~,Jg]=hexa(A,Xloc,u0loc); % create local matrix and rhs
  61            %disp('Jg is')
  62            %disp(Jg)
  63             % instead, accumulate contributions to the global matrix
  64             if ~isempty(lambda)
  65                 grad = lambda(kglo)'*Jg;  % grad lambda
  66                %disp('Local lambda is')
  67                %lambda(kglo)
  68                 %disp(' Grad prior to multiplication by A_inv is')
  69                %disp(grad)
  70                 grad = grad/A;
  71                 %disp('Grad after to multiplication by A_inv is')
  72                  %disp(grad)
  73                 for i=1:3
  74                     W{i}(i1,i2,i3)=grad(i);
  75                 end
  76             end
  77         end
  78     end
  79 end
  80 if ~isempty(u0)
  81     for i=1:3
  82         W{i}=u0{i}+W{i};
  83     end
  84 end
  85 % fprintf('\n')
  86 end