% run capacitormg
clear all
global Prmat method presm postsm relax lmax
format compact
% size of inner conductor
a   = .75; 
b   = 0.5; 
% parameters: number of gridpoints 
m2  = 9;
n   = 2^m2;
tol = 1e-8;
% what is the optimal relaxtion factor?
relax = .5;
% number of MG levels (level 0 = finest)
lmax = log2(n)-3,
method = 'MG2';  % Jacobi or MG2 
presm  = 2;
postsm = 2;
tic  ;           % Start timer
[err,capT,WU,pot,mask]=capacitormg(a,b,n,tol,relax);
csp = -(log10(err(end))-log10(err(1)))/WU;
nit=length(err);
fprintf(['n: %3i, Method: %s, # iterations:%4i, #digits/WU:' ...
					   ' %5.3f \n',n,method,nit,csp);
toc;                        % stop timer
subplot(221); semilogy(err);title(['//fn+1-fn//2, n = ', ...
		    num2str(n)])
subplot(222); mesh(pot); title('V(x,y)');
subplot(223); semilogy(abs(capT(1:end-1)-capT(end))), ...
    title('cap(n+1)-cap(n)');

------------------------------

function [err,capT,WU,f,mask] = capacitormg(a, b, n, tol, rel)
% --------------------------------------------------------------
% Compute capacitance per unit length of
% a coaxial pair of rectangles
% --------------------------------------------------------------
global Prmat method presm postsm relax mask lmax
% Arguments:
%    a   =  width of inner conductor
%    b   = height of inner conductor
%    n   = number of points in the x-direction (horizontal)
%    tol = relative tolerance for capacitance
%    rel = relaxation parameter

% Returns:
%    capT(1..nit) = capacitance per unit length 
%    err(1..nit)    list of ||fn+1 - fn||2
%    WU             total # work units
%    f(1:n+1,1:n+1) potential
%    mask(d:o)      = 1 between conductors
%
clear capT err gerr

% Make grids: enclosing unit square
x  = linspace(-0.5,0.5,n+1);  % Grid points along x-axis
m  = n;                  % Number of segments on 'd'
y  = x;  % Grid points along y-axis
[xx,yy]=meshgrid(x,y);       % 2D grid

% set up prolongation matrices
mm    = max(m+1,n+1);
% Prmat = ...

% Initialize potential and mask array
% [mask(i,j) = 0 implies unchanged f(i,j)]
f    = zeros(n+1,m+1);       % 2D-array with solution
mask = ones(n+1,m+1)*rel;    % 2D-array with relaxation factor

% enclosing PEC rectangle:
mask(1:n+1,m+1)=0;
mask(1:n+1,  1)=0;
mask(n+1,1:m+1)=0;
mask(  1,1:m+1)=0;

% inner conductor
geom = 'square';
geom = 'ellipse';
switch geom
    case 'square'
        inside = find((abs(xx)<=(a/2))&(abs(yy)<=(b/2)));
    case 'ellipse'
        inside = find((xx/(a/2)).^2 + (yy/(b/2)).^2 <= 1);
end
mask(inside)=0;
f(inside)=1;

oldcap = -10;
WU = 0;
maxit = 20;
for iter = 1:maxit         % Maximum number of iterations
    f0 = f;
    switch method          % Perform iteration
         case 'Jacobi'
            f = jacobi2(f,0*ones(size(f))/n^2,mask);
            WU = WU+1;
        case 'MG2'
            f  = MGit2(f,0*ones(size(f))/n^2,mask,0,lmax);
            WU = WU+(presm+postsm+1)*(1+1/3);
    end
    err(iter)=norm(f(:)-f0(:),2);
    disp(err(iter));
    cap = gauss(f);               % Compute the capacitance
    capT(iter) = cap;
    if (abs(cap-oldcap)<tol*cap)
        break                     % Stop if change in capacitance
                                  % is sufficiently small
    else
        oldcap = cap;             % Continue
    end
end

------------------------------
function u = MGit2(u,rhs,mask,l,lmax);
global Prmat presm postsm relax
[n,m]=size(u);

for k = 1:presm
    u = jacobi2(u,rhs,mask);
end
if l < lmax
% compute prolongation and restriction matrices
%    Pl, Pr, Rl, Rr
% compute residual
    res = (jacobi2(u,rhs,mask)-u)/relax;
%    subplot(4,lmax,l+1);
%    mesh(res)
%    title(num2str(l))
% solve correction equation by sending restricted residual to
% coarser grid:
    u2 = MGit2(zeros(size(u(1:2:end,1:2:end))),Rl*res*Rr,mask(1:2: ...
						  end,1:2:end),l+1,lmax);
% prolong correction
    corr = Pl*u2*Pr;
%    subplot(4,lmax,l+1+lmax); mesh(corr);  
% add
    u  = u+corr;
%    subplot(4,lmax,l+1+2*lmax); mesh(u);
    % reset Dirichlet conditions
    inside    = find(mask==0);
    u(inside) = (l==0);
    u(n,:)=0;
    u(1,:)=0;
    u(:,m)=0;
    u(:,1)=0;
end
kk = postsm;
if l==lmax
% .. solve more exactly on coarsest grid
    kk = 10;
end
for k = 1:kk
    u = jacobi2(u,rhs,mask);
end
if l < lmax 
    %subplot(4,lmax,l+1+3*lmax); mesh(u);pause 
end

--------------------------------------
function u = jacobi2(u,rhs,mask)
% Arguments:
%    u    = 2D-array with solution
%    rhs  = 2D array with right hand side, res = laplace u + rhs
%    mask = 2D-array with relaxation/time step factor 
% Returns:
%    u    = 2D-array with solution after Jacobi iteration, i.e.,
%  u(n+1)ij = u(n)ij + relax*maskij*(A*u(n)+rhs)ij
% Tips: A = d2/dx2 + d2/dy2; expand u by rows & columns to make
% diff(uexp,2,k), k = 1 and 2, work;
% The ulterior rows & columns can be anything because those results
% are masked out 
% compute "Au := A*u"
% update
u   = u+mask.*(Au+rhs);

---------------------------------------------
function q = gauss(f)
% --------------------------------------------------------------
% Compute capacitance from the potential
% Assumes dx = dy!
% --------------------------------------------------------------
% Arguments:
%    n    = number of points in the x-direction (horizontal)
%    m    = number of points in the y-direction (vertical)
%    f    = 2D-array with solution
% Returns:
%    q     = capacitance per unit length 
% uses curve close to outer boundary
% modify to other curve, checlk that this gives same result
[n1,m1] = size(f);
n = n1-1;
m = m1-1;
q = 0;

for i = 1:n
  q = q + (f(i,m)+f(i+1,m))*0.5+...
          (f(i,2)+f(i+1,2))*0.5 ; %integrate along
                                  %horiz. boundaries
end
for j = 1:m
  q = q + (f(n,j)+f(n,j+1))*0.5+...
          (f(2,j)+f(2,j+1))*0.5; %integrate along
                                 %vertic. boundaries
end




