leapfrog.m

解决麦克斯韦的matlab源文件

function [energies,sol_v,sol_e,times] = leapfrog(mesh,vectfld,ts,T,grabstep,scal)
% leapfrog timestepping for two discretizations of Maxwell's equations
%
% mesh -> 2D unstructured mesh
% vectfld -> string designating the routine providing the initial vectorfield
% rs -> timestep
% T -> final time
% grabstep -> every #grabstep iterate will be sampled
% scal ->  governs strength of regularization A = 2*(scal*C + (1-scal)*R)
%
% Result:
% 
% energies -> trace of electric/magnetic energies during timestepping
%        energies(:,1) = time, 
%        energies(:,2) = electric energy of nodal solution
%        energies(:,3) = magnetic energy of nodal solution
%        energies(:,4) = electric energy of edge element solution
%        energies(:,5) = magnetic energy of edge element solution
% sol_v -> sampled solutions for nodal FEM
% sol_e -> sampled solutions for edge elements
% times -> vector of sampling times
%

if (nargin < 6), scal = 0.5; end
if (nargin < 5), grabstep = 1; end

% Here the crucial matrices are obtained. 
% NOTE: Different scalings of curl-curl and regularizing part
% should be taken into account here.
% See documentation of Emats.m/Nmats.m

disp('Assembling matrices');
[Ae,Me] = Emats(mesh,1,scal);
[An,Mn] = Nmats(mesh,1,scal);

% Initial values for nodal scheme
% Plain nodal interpolation for vertex based elements

disp('Computing initial field for nodal scheme');
v_vals = p1_interpolate(vectfld,mesh);
n_act = Nactidx(mesh);
nv = v_vals(n_act);

% Initial values for edge element scheme
% Interpolation + projection

disp('Computing E_0 for edge elements');
e_vals = ed_interpolate(vectfld,mesh);
e_act = Eactidx(mesh);
ev = e_vals(e_act);
G_all = Gradmat(mesh);
v_act = Vactidx(mesh);
G = G_all(e_act,v_act);
p = G*((G'*Me*G)\(G'*Me*ev));
fprintf('Norm of edge interpolant = %f, norm of correction = %f\n',...
       norm(ev),norm(p));
ev = ev - p;
clear e_vals v_vals n_act e_act G_all G p

nv_new = zeros(size(nv));
nv_mid = zeros(size(nv));
nv_tmp = zeros(size(nv));
ev_new = zeros(size(ev));
ev_mid = zeros(size(ev));
ev_tmp = zeros(size(ev));

disp('Displaying initial iterates');
figure; %figure('Position',get(0,'ScreenSize')); 
clf; figno = gcf;
plotiterate(mesh,ev,nv,0,figno);

sol_v = nv;
sol_e = ev;
times = [0.0];

%F(1) = getframe(figno);

nv_tmp = An*nv;
ev_tmp = Ae*ev;
etot_v = dot(nv_tmp,nv);
etot_e = dot(ev_tmp,ev);

disp('Initial energies:');
fprintf('\t #### Nodal scheme : E_el = %f, E_mag = %f\n',...
	0,etot_v);
fprintf('\t #### Edge elements: E_el = %f, E_mag = %f\n',...
	0,etot_e);

% First step
nv_mid = nv - 0.5*ts*ts*(Mn\nv_tmp);
ev_mid = ev - 0.5*ts*ts*(Me\ev_tmp);

stp = 1;
if (grabstep == 1)
  sol_v = [sol_v nv_mid];
  sol_e = [sol_e ev_mid];
  times = [times ts];
end

plotiterate(mesh,ev_mid,nv_mid,ts,gcf);

energies = [0.0 0.0 etot_v 0.0 etot_e];
for t=2*ts:ts:T
  stp = stp + 1;
  fprintf('Iteration step %d at time %d\n',stp,t);
  nv_tmp = An*nv_mid;
  ev_tmp = Ae*ev_mid;
  nv_new = 2*nv_mid - nv - ts*ts*(Mn\nv_tmp);
  ev_new = 2*ev_mid - ev - ts*ts*(Me\ev_tmp);
  disp('Displaying new iterate');
  plotiterate(mesh,ev_new,nv_new,t,figno);
  if (mod(stp,grabstep) == 0)
    sol_v = [sol_v nv_new];
    sol_e = [sol_e ev_new];
    times = [times t];
    %F(stp/grabstep) = getframe(figno); 
  end
  
  % Computing energies
  nv = (nv_new - nv)/(2*ts);
  ev = (ev_new - ev)/(2*ts);
  el_en_v = dot(nv,Mn*nv);
  el_en_e = dot(ev,Me*ev);
  mag_en_v = dot(nv_tmp,nv_mid);
  mag_en_e = dot(ev_tmp,ev_mid);
  
  fprintf('\t Nodal scheme : E_el = %f, E_mag = %f, E_tot = %f\n',...
	  el_en_v,mag_en_v,el_en_v+mag_en_v);
  fprintf('\t Edge elements: E_el = %f, E_mag = %f, E_tot = %f\n',...
	  el_en_e,mag_en_e,el_en_e+mag_en_e);
  if (el_en_v+mag_en_v > 10*etot_v)
    disp('Instability of nodal scheme!');
    break;
  end
  if (el_en_e+mag_en_e > 10*etot_e)
    disp('Instability of edge element scheme!');
    break;
  end
  
  energies = [energies; t el_en_v mag_en_v el_en_e mag_en_e];
  nv = nv_mid;
  ev = ev_mid;
  ev_mid = ev_new;
  nv_mid = nv_new;
end