nmats.m

解决麦克斯韦的matlab源文件

function [A,M] = Nmats(mesh,withbd,scal)
% computes the matrices A ("stiffness matrix") and M ("mass matrix")
% for continuous piecewise linear vectorial finite elements
%
% Usage: [A,M] = Nmats(mesh)
%
% mesh -> data structure for 2D triangulation
% withbd -> If true homogeneous b.c. are assumed, default TRUE
% scal -> strength of regularization: A = 2*scal*C + 2*(1-scal)*R
%

if (nargin < 2), withbd = 1; end
if (nargin < 3), scal = 0.5; end

disp('Assembling nodal element matrices');
fprintf('Relative scaling of C and R: %d <-> %d\n',2*scal,2*(1-scal));

% Dimension of space without boundary conditions
N = 2*mesh.Nv;

Ap = sparse(N,N);
Mp = sparse(N,N);

% Local assembly procedure: Run through all triangles
% and sum up element matrices

for idx=mesh.trv'
  x = mesh.vt(idx,1);
  y = mesh.vt(idx,2);
  
% Compute the determinant and the area.
  d = (x(2)-x(1))*(y(3)-y(1))-(x(3)-x(1))*(y(2)-y(1));
  area = abs(d)/2;
  
% Local mass matrix
  MT = area*[2,0,1,0,1,0;...
	     0,2,0,1,0,1;...
	     1,0,2,0,1,0;...
	     0,1,0,2,0,1;...
	     1,0,1,0,2,0;...
	     0,1,0,1,0,2]/12.0;
% Local curl-curl matrix
  K = [ x(3)-x(2),y(3)-y(2),x(1)-x(3),y(1)-y(3),x(2)-x(1),y(2)-y(1) ]/(2*area);
  CT = area*K'*K;
  
% Local grad-div matrix
  L = [ y(2)-y(3),x(3)-x(2),y(3)-y(1),x(1)-x(3),y(1)-y(2),x(2)-x(1) ]/(2*area);
  RT = area*L'*L;

% Local stiffness matrix
  AT = 2*scal*CT + 2*(1-scal)*RT;
  
% Global numbers of degrees of freedom
% Basis functions for x-components => odd indices
% Basis functions for y-components => even indices
  dofs = [2*idx(1)-1,2*idx(1),2*idx(2)-1,2*idx(2),2*idx(3)-1,2*idx(3)];
  
  Ap(dofs,dofs) = Ap(dofs,dofs) + AT;
  Mp(dofs,dofs) = Mp(dofs,dofs) + MT;
end

% Elimination of contributions of inactive vertices on the boundary

act_idx = Nactidx(mesh);

if (withbd == 1)
  A = Ap(act_idx,act_idx);
  M = Mp(act_idx,act_idx);
else
  A = Ap;
  M = Mp;
end