gm_assemble.m

matlab有限元分析工具,比经较全面的一个手册,请大家下载呀

function [K,f1] = gm_assemble(vtx, tri, isdir, dirval, conduc)
% [K,f1] = gm_assemble(vtx, tri, isdir, dirval, conduc)
% This function assembles the finite element stiffness matrix K
% and the portion of the rhs f that is due to Dirichlet BC's.
%  vtx: an n-by-d zba matrix of vertex coordinates
%  tri: an m-by-(d+1) zba of simplex vertices (integer indices into vtx).
%  isdir: tells which vertices have dirichlet BC's.  
%       An n-vector (zba) of 1's and 0's.
%  dirval: An n-vector (zba) of dirichlet values.  Only positions in which
%   isdir = 1 are used.
%  conduc: Conductivity field (m-by-1 array of positive reals) averaged
%    over each simplex, zba.
%  This function is called by gmfem.

[numnodes, di] = size(vtx);
[numtri, scrap] = size(tri);

% The goal of this first section is to make a matrix A such that
% if A is multiplied by a vector of nodal values (all nodes except
% dir. BC's) then it returns the element-by-element gradient (with
% di components per gradient).  Thus, A is (m*di)-by-n1 where m is
% the number of elements, di is the dimension, and n1 is the number
% of nonboundary nodes.

% Make a combinatorial version of A that has 1's, 0's and -1's.
% Call it A0.

nondir = zba(ones(numnodes,1)) - isdir;

ilist = kron((1:di*numtri), [1,1]);
jlist = zeros(1,numtri*di*2);
slist = zeros(1,numtri*di*2);

for d = 1 : di
  jlist(2 * d - 1 : 2 * di : numtri * di * 2) = double(tri(:,0))'+...
          ones(1,numtri);
  jlist(2 * d : 2 * di : numtri * di * 2) = double(tri(:,d))' + ...
          ones(1,numtri);
  slist(2 * d - 1 : 2 * di : numtri * di * 2) = -ones(1,numtri);
  slist(2 * d : 2 * di : numtri * di * 2) = ones(1,numtri);
end

A0 = sparse(ilist,jlist,slist,di*numtri,numnodes);


%%%% Make A0a be A0, with columns corresponding to
%%%% nonboundaries deleted.

findisdir = find(isdir);
findisdir1 = double(findisdir) + ones(length(findisdir),1);

A0a = A0(:,findisdir1);


%%%% Make A0b be A0, with columns corresponding to
%%%% boundaries deleted.  A will be based on A0b.

findnondir = find(nondir);
findnondir1 = double(findnondir) + ones(length(findnondir),1);
A0b = A0(:,findnondir1);

% Make the block diagonal matrix iJ with the Jacobians of the
% elements.  iJ is composed of m blocks, each block being di-by-di.

if di == 2
  J1 = A0 * double(vtx(:,0));
  J2 = A0 * double(vtx(:,1));
  J11 = J1(1:2:numtri*2);
  J12 = J1(2:2:numtri*2);
  J21 = J2(1:2:numtri*2);
  J22 = J2(2:2:numtri*2);
  detJ = J11 .* J22 - J21 .* J12;
  % msl = sqrt(max([(J11.^2+J21.^2),(J12.^2+J22.^2)]'))';
  invJ11 = J22 ./ detJ;
  invJ12 = -J12 ./ detJ;
  invJ21 = -J21 ./ detJ;
  invJ22 = J11 ./ detJ;
elseif di == 3
  J1 = A0 * double(vtx(:,0));
  J2 = A0 * double(vtx(:,1));
  J3 = A0 * double(vtx(:,2));
  J11 = J1(1:3:numtri*3);
  J12 = J1(2:3:numtri*3);
  J13 = J1(3:3:numtri*3);
  J21 = J2(1:3:numtri*3);
  J22 = J2(2:3:numtri*3);
  J23 = J2(3:3:numtri*3);
  J31 = J3(1:3:numtri*3);
  J32 = J3(2:3:numtri*3);
  J33 = J3(3:3:numtri*3);
  detJ = J11 .* J22 .* J33 - J11 .* J23 .* J32 - J12 .* J21 .* J33 ...
          + J12 .* J23 .* J31 + J13 .* J21 .* J32 - J13 .* J22 .* J31;
  %msl = sqrt(max([(J11.^2+J21.^2+J31.^2),(J12.^2+J22.^2+J32.^2), ...
	% (J13.^2+J23.^2+J33.^2)]'))';
  invJ11 = (J22 .* J33 - J23 .* J32) ./ detJ;
  invJ12 = (J32 .* J13 - J12 .* J33) ./ detJ;
  invJ13 = (J12 .* J23 - J22 .* J13) ./ detJ;
  invJ21 = (J31 .* J23 - J21 .* J33) ./ detJ;
  invJ22 = (J11 .* J33 - J13 .* J31) ./ detJ;
  invJ23 = (J21 .* J13 - J11 .* J23) ./ detJ;
  invJ31 = (J21 .* J32 - J31 .* J22) ./ detJ;
  invJ32 = (J31 .* J12 - J11 .* J32) ./ detJ;
  invJ33 = (J11 .* J22 - J21 .* J12) ./ detJ;
else
  error(' Only dim 2 and 3 supported.')
end


ilist = kron((1 : di * numtri), ones(1,di));
jlist = zeros(1, numtri*di^2);
for d = 1 : di 
  jlist(d:di:numtri*di^2) = kron((d:di:di*numtri),ones(1,di));
end

if di == 2
  slist = zeros(1,numtri*4);
  slist(1:4:numtri*4) = invJ11;
  slist(2:4:numtri*4) = invJ21;
  slist(3:4:numtri*4) = invJ12;
  slist(4:4:numtri*4) = invJ22;
else
  slist = zeros(1,numtri*9);
  slist(1:9:numtri*9) = invJ11;
  slist(2:9:numtri*9) = invJ21;
  slist(3:9:numtri*9) = invJ31;
  slist(4:9:numtri*9) = invJ12;
  slist(5:9:numtri*9) = invJ22;
  slist(6:9:numtri*9) = invJ32;
  slist(7:9:numtri*9) = invJ13;
  slist(8:9:numtri*9) = invJ23;
  slist(9:9:numtri*9) = invJ33;
end


  % slist = slist .* kron(msl',[1,1,1,1]);
iJ = sparse(ilist,jlist,slist,di*numtri,di*numtri);


% A is the "reduced geometry matrix" described above.

A = iJ * A0b;

% get the dirichlet bc's in b

b = - iJ * (A0a * double(dirval(find(isdir))));

% get the areas of the elements, which is |detJ|/di!.

vols = abs(detJ) / prod(1:di);

% Get the diagonal weight matrix D that has conductivity integrated.

iD = sparse((1:di*numtri), (1:di*numtri),  ...
    kron((double(conduc).*vols)', ones(1,di)));

K = A' * iD * A;
f1 = A' * iD * b;


% ------------------------------------------------------------------
% Copyright (c) 1998 by Cornell University.  All rights reserved.
% See the accompanying file 'Copyright' for authorship information,
% the terms of the license governing this software, and disclaimers
% concerning this software.
% ------------------------------------------------------------------
% This file is part of the QMG software.  
% Version 2.0 of QMG, release date RELDATE
% ------------------------------------------------------------------