femwind/nd_assembly.m

   1 function K=nd_assembly(A,X,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 % out:
   6 %  K  stiffness matrix, stored as size (n1,n2,n3,3,3,3)
   7
   8 fprintf('sparse_assembly:')
   9
  10 d = size(X,2);    % dimensions
  11 n = size(X{1});   % mesh size in nodes
  12 nn = prod(n);     % total nodes
  13
  14 % initialize matrices
  15 F = zeros(nn,1);
  16 K = zeros(n(1),n(2),n(3),3,3,3);
  17
  18 % preallocate matrices
  19 tstart=tic;    % start timer
  20 Xloc = zeros(3,8);
  21 kglo=zeros(1,8);
  22 m = n-1;       % grid size in elements
  23 disp('accumulating global stiffness matrix entries')
  24 for i3=1:m(3)
  25     for i2=1:m(2)
  26         for i1=1:m(1)  % loop over elements
  27             % now in element (i1,i2,i3)
  28             % build the matrix of node coordinates to pass to hexa
  29             for j3=0:1 % loop over corners of the element
  30                 for j2=0:1
  31                     for j1=0:1
  32                         jloc=1+j1+2*(j2+2*j3);  % local index of the node in element (i1, i2. i3)
  33                         for i=1:3
  34                             Xloc(i,jloc)=X{i}(i1+j1,i2+j2,i3+j3); % node coordinate i
  35                         end
  36                     end
  37                 end
  38             end
  39             [Kloc,~,~]=hexa(A,Xloc,zeros(3,1)); % compute the local matrix
  40             % loop over local matrix entries j,k
  41             for j3=0:1 %
  42                 for j2=0:1
  43                     for j1=0:1
  44                         jloc=1+j1+2*(j2+2*j3); % local index j
  45                         for k3=0:1 % loop over corners of the element
  46                             for k2=0:1
  47                                 for k1=0:1
  48                                     kloc=1+k1+2*(k2+2*k3);
  49                                     % add entry of the local matrix; need
  50                                     % to offset by 2 because all arrays
  51                                     % dimensions in matlab start at 1
  52                                     K(i1+j1,i2+j2,i3+j3,2+k1-j1,2+k2-j2,2+k3-j3)=...
  53                                         K(i1+j1,i2+j2,i3+j3,2+k1-j1,2+k2-j2,2+k3-j3)+...
  54                                         Kloc(jloc,kloc);
  55                                 end
  56                             end
  57                         end
  58                     end
  59                 end
  60             end
  61         end
  62     end
  63 end
  64
  65 end