Skip to content

Instantly share code, notes, and snippets.

@ashgillman
Created October 25, 2015 10:11
Show Gist options
  • Save ashgillman/58c6273628d4e73301af to your computer and use it in GitHub Desktop.
Save ashgillman/58c6273628d4e73301af to your computer and use it in GitHub Desktop.
Laplace Finite Differences
clear all; close all
% Input Parameters
shieldxs = [0, 10, 10*cos(pi/3), 0];
shieldys = [0, 0, 10*sin(pi/3), 0];
shieldv = 10;
corexs = [4, 4, 6, 6, 4];
coreys = [5*tan(pi/6)+1, 5*tan(pi/6)-1, 5*tan(pi/6)-1, ...
5*tan(pi/6)+1, 5*tan(pi/6)+1];
corev = -10;
%gridsize = [129, 129];
maxits = 10e6;
err = 10e-6;
exponentRange = 3:15;
sampleV = [0.5*max(shieldxs), 0.75*max(shieldys)]; %sample voltage in real units
Vs = zeros(1, length(exponentRange));
for ei=1:length(exponentRange)
e = exponentRange(ei);
gridsize = [2^e+1, 2^e+1];
% Field (real units)
realsize = [max(shieldxs), max(shieldys)];
xs = linspace(0, realsize(1), gridsize(1));
ys = linspace(0, realsize(2), gridsize(2));
dx = realsize(1)/(gridsize(1)-1);
dy = realsize(2)/(gridsize(2)-1);
units2grid = @(x, d) x./realsize(d).*(gridsize(d)-1);
% Grid (pixels)
Vgrid = ones(gridsize);
gridxs = repmat(0:gridsize(1)-1, 1, gridsize(2));
gridys = repmat(0:gridsize(2)-1, gridsize(1), 1);
gridys=gridys(:)';
% Set grid initial values
Vgrid(:) = (shieldv+corev) / 2; % set initial est V
[in,on] = inpolygon(gridxs', gridys', ...
units2grid(shieldxs, 1), units2grid(shieldys, 2));
shieldpts = ~in|on;
Vgrid(shieldpts) = shieldv; % set shield V
[in,on] = inpolygon(gridxs, gridys, ...
units2grid(corexs, 1), units2grid(coreys, 2));
corepts = in;
Vgrid(corepts) = corev; % set core V
% iterations
i = 2:gridsize(1)-1;
j = 2:gridsize(2)-1;
for it=1:maxits
old = Vgrid;
Vgrid(i,j) = 0.5*(dx*(Vgrid(i, j-1)+Vgrid(i, j+1)) + ...
dy*(Vgrid(i-1, j)+Vgrid(i+1, j))) / (dx+dy);
Vgrid(shieldpts) = shieldv; % reset shield V
Vgrid(corepts) = corev; % reset core V
if (max(max(abs(Vgrid-old))) < err) % convergence
break
end
end
disp(['broke after ' num2str(it) ' its']);
Vs(ei) = Vgrid(units2grid(sampleV(1), 1), units2grid(sampleV(2), 2));
disp(['V at (' num2str(sampleV(1)) ', ' num2str(sampleV(2)) ...
') is ' num2str(Vs(ei))]);
% Calculate E
Ex = zeros(size(Vgrid));
Ey = zeros(size(Vgrid));
Ex(i,j) = (Vgrid(i-1,j)-Vgrid(i+1,j)) / (2*dx);
Ey(i,j) = (Vgrid(i,j-1)-Vgrid(i,j+1)) / (2*dy);
Ex(shieldpts) = 0; % reset shield E
Ey(shieldpts) = 0;
Ex(corepts) = 0; % reset core E
Ey(corepts) = 0;
E = sqrt(Ex.^2 + Ey.^2);
% Plot V and E
down = 4;
Ex = downsample(downsample(Ex, down)', down)'; % Downsample E for plot
Ey = downsample(downsample(Ey, down)', down)';
figure
hold on
surf(xs, ys, Vgrid','LineStyle','none');
shading interp; colorbar
line(shieldxs,shieldys,max(max(Vgrid))*ones(size(shieldxs)), ...
'color','w');
line(corexs,coreys,max(max(Vgrid))*ones(size(corexs)),'color','w');
axis tight; axis equal
title('V field'); xlabel('x (cm)'); ylabel('y (cm)')
figure
hold on
quiver(downsample(xs, down), downsample(ys, down), Ex', Ey');
contour(xs, ys, E', 4);
colorbar
plot(shieldxs, shieldys, 'color','k');
plot(corexs, coreys, 'color','k');
axis tight; axis equal
title('E field: contours show lines of constant |E|');
xlabel('x (cm)'); ylabel('y (cm)')
figure(3)
plot(10./2.^exponentRange(1:ei), Vs(1:ei), 'o-')
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment