mptransfer.m

解决麦克斯韦的matlab源文件

function mpv = mptransfer(mesh,nv)
% Computes values of edge interpolant at midpoints
%
% mesh -> data structure for 2D unstructured grid
% nv -> column vector of legth 2*mesh.Nv, containing x/y-components of 
%       vectorfield at vertices
%
% Result ev -> column vector containing values at midpoints 
%

if (size(nv,1) ~= 2*mesh.Nv), error('Size mismatch for argument vector'); end
if (size(nv,2) ~= 1), error('nv must be a column vector'); end

mpv = zeros(mesh.Nt,2);

for i=1:mesh.Nt
  vidx = mesh.trv(i,:);
  x = mesh.vt(vidx,1);
  y = mesh.vt(vidx,2);
  e1 = [x(3)-x(2); y(3)-y(2)];
  e2 = [x(1)-x(3); y(1)-y(3)];
  e3 = [x(2)-x(1); y(2)-y(1)];
  v1 = [nv(2*vidx(1)-1); nv(2*vidx(1))];
  v2 = [nv(2*vidx(2)-1); nv(2*vidx(2))];
  v3 = [nv(2*vidx(3)-1); nv(2*vidx(3))];
  dof1 = 0.5*dot(v2+v3,e1);
  dof2 = 0.5*dot(v1+v3,e2);
  dof3 = 0.5*dot(v2+v1,e3);
  
  disp('Geometry'); disp([x';y']);
  disp('Edges'); disp([e1 e2 e3]);
  disp('Vertex values'); disp([v1 v2 v3]);
  disp('DOF'); disp([dof1 dof2 dof3]);
  d = (x(2)-x(1))*(y(3)-y(1))-(x(3)-x(1))*(y(2)-y(1));
  area = abs(d)/2;
  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)];
  val = [G(:,3)-G(:,2) , G(:,1)-G(:,3) , G(:,2)-G(:,1)]*[dof1;dof2;dof3]/3;
  disp('Gradients'); disp(G);
  mpv(i,:) = val';
  disp('Central value'); disp(val);
  pause;
end