emats.m

解决麦克斯韦的matlab源文件

function [A,M] = Emats(mesh,withbd,scal)
% Computes "mass matrix" and (discretely regularized) "stiffness matrix"
% for edge elements
%
% mesh -> data structure for 2D triangulation
% withbd -> If true drop d.o.f. on boundary, default = TRUE
% scal -> governs strength of regularization A = 2*(scal*C + (1-scal)*R)
%


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

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

N = mesh.Ne;

Ap = sparse(N,N);
Mp = sparse(N,N);
Dp = zeros(mesh.Nv,1);

% Assembly of local matrices

for i=1:mesh.Nt
  vidx = mesh.trv(i,:);
  x = mesh.vt(vidx,1);
  y = mesh.vt(vidx,2);
  dofs = mesh.tre(i,:);
  
% 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;

% Gradients of barycentric coordinate functions
  G = 1/(2*area)*[y(2)-y(3) , y(3)-y(1) , y(1)-y(2);...
		  x(3)-x(2) , x(1)-x(3) , x(2)-x(1)];
  
% Local mass matrix
  P = 0.5*[ G(:,3)-G(:,2) ,   G(:,1)      ,  -G(:,1)     ; ...
	    -G(:,2)       , G(:,1)-G(:,3) ,  G(:,2)      ;
	    G(:,3)        ,    -G(:,3)    , G(:,2)-G(:,1) ];
  MT = area/3*P'*P;
  
% Local curl-curl matrix
  Q = [1,1,1]/area;
  CT = area*Q'*Q;

% Get permutation matrix reflecting the mismatch
% between global and local edge orientations
  Peori = diag(mesh.treo(i,:));
  
  Ap(dofs,dofs) = Ap(dofs,dofs) + Peori*CT*Peori;
  Mp(dofs,dofs) = Mp(dofs,dofs) + Peori*MT*Peori;
  
% Diagonal lumped mass matrix for p.w. linear Lagrangian
% finite elements

  DT = area/3*[1;1;1];
  Dp(vidx',1) = Dp(vidx',1) + DT;
end

% Topological gradient matrix
G = Gradmat(mesh);

% Add the two contributions to total stiffness matrix

Am = 2*scal*Ap + 2*(1-scal)*G*spdiags(Dp,0,mesh.Nv,mesh.Nv)*G';

% Discard rows and columns corresponding to edges on the boundary

if (withbd == 1)
  act_idx = Eactidx(mesh);
  A = Am(act_idx,act_idx);
  M = Mp(act_idx,act_idx);
else
  A = Am;
  M = Mp;
end;