femwind/apply_boundary_conditions.m

   1 function [Kb,Fb]=apply_boundary_conditions(K,F,X)
   2 % apply zero dirichet boundary condition on sides and on top
   3 % in:
   4 %   K   stiffness matrix
   5 %   F   load vector
   6 %   X   X{1} has the size of the domain
   7 % out:
   8 %   Kb   updated
   9 %   Fb   updated
  10
  11 n = size(X{1});
  12 nn=prod(n);
  13 [i1,i2,i3]=ind2sub(n,1:nn);
  14 bc = i1==1 | i1 ==n(1) | i2 ==1 | i2 == n(2) | i3 == n(3);  % sides and top only, bottom left alone
  15 bx = sub2ind(n,i1(bc),i2(bc),i3(bc)); % matrix indices with zero boundary condition
  16 % check number of nodes with boundary condition
  17 nbc = (n(1)-2)*(n(2)-2)+2*(n(1)-1)*n(3)+2*(n(2)-1)*n(3);
  18 if nnz(bc) ~= nbc,
  19     error('wrong number of boundary nodes')
  20 end
  21 Kb=[]; Fb=[];
  22 if ~isempty(K)
  23     if any(size(K)~=nn)
  24         error('apply_boundary_conditions: inconsistent size K')
  25     end
  26     Kb=K;
  27     Kb(bx,:)=0;
  28     Kb(:,bx)=0;
  29     Kb(bx,:)=0;
  30     Kb(bx,bx)=big(diag(K))*speye(length(bx));
  31 end
  32 if ~isempty(F)
  33     Fb=F;
  34     Fb(bx)=0;
  35 end
  36 end