femwind/f_assembly.m

   1 function F=f_assembly(A,X,u0,params);
   2
   3 % d = size(X,2);    % dimensions
   4 n = size(X{1});   % mesh size in nodes
   5 nn = prod(n);     % total nodes
   6
   7
   8 % initialize matrices
   9 F = zeros(n(1),n(2),n(3));
  10 Xloc = zeros(3,8);
  11 kglo=zeros(1,8);
  12 u0loc = zeros(3);
  13
  14 for ie2=1:n(2)-1  % loop over elements
  15     for ie3=1:n(3)-1
  16         for ie1=1:n(1)-1
  17             % now in element (ie1,ie2,ie3)
  18             % build the matrix of node coordinates to pass to hexa
  19             for ic2=0:1 % loop over corners of the element
  20                 for ic3=0:1
  21                     for ic1=0:1
  22                         iloc=1+ic1+2*(ic2+2*ic3);  % local index of the node in the element
  23                          %  position of the node in the global grid
  24                         %kglo(iloc)=k1+n(1)*((k2-1)+n(2)*(k3-1)); % global index
  25                         for i=1:3
  26                             Xloc(i,iloc)=X{i}(ie1+ic1,ie2+ic2,ie3+ic3); % node coordinate i
  27                         end
  28                     end
  29                 end
  30             end
  31             if ~isempty(u0)
  32                 u0loc =[u0{1}(ie1,ie2,ie3),u0{2}(ie1,ie2,ie3),u0{3}(ie1,ie2,ie3)]';
  33             else
  34                 u0loc=[];
  35             end
  36
  37             [~,Floc,~]=hexa(A,Xloc,u0loc); % compute the local load vector
  38             for id2 = 0:1
  39                 for id3 = 0:1
  40                     for id1 = 0:1
  41                         iloc=1+id1+2*(id2+2*id3);
  42                         k1 = ie1+id1; k2 = ie2+id2; k3 = ie3+id3;
  43                         F(k1,k2,k3) = F(k1,k2,k3) + Floc(iloc);
  44                     end
  45                 end
  46             end
  47         end
  48     end
  49 end
  50 end