% The 8 fields "XXX" below should be specified

%
% This is a Matlab script FEM.m which produces a 3d plot of the solution to the Dirichlet problem, which it calculates using the finite element method which mesh size 1/L and using bilinear elements on squares. It requires Matlab programs h1.m and h2.m for the lower and upper boundaries of the domain and f.m for (an extension of) the boundary data
%
%

L= XXX ; % mesh size =1/L
a= XXX ; % x-coordinate for left end of domain 
b= XXX ; % x-coordinate for right end of domain 
% Required functions: 
%   lower boundary curve y=h1(x)
%   upper boundary curve y=h2(x)
%   extension of boundary Dirichlet data z=f(x,y)

dx= (b-a)/(L+1); 
y1= zeros(L,1);
y2=zeros(L,1);
for k=1:L
  y1(k)=floor(h1(a+k*dx)/dx);  %lower boundary y1= integer valued scale
  y2(k)=floor(h2(a+k*dx)/dx);  %upper boundary y2= integer valued scale
end
ydiff= y2-y1+1; %y1= integer valued length of vertical cross section

N=sum(ydiff);  % total number of trial functions

M=zeros(N,N);  % stiffness matrix M
G=zeros(N,1);  % source vector G
F=zeros(N,1);  % extension of Dirichlet data vector F

i1=0;   % start index of current column in the domain
for n=1:L  % column in the domain
    i2=i1+ydiff(n); % start index of next column in the domain 
    for j=1:ydiff(n)
        % neighbours in the stiffness matrix
        if j<ydiff(n)
            M(i1+j,i1+j+1)=-1; %upper neighbour
        end
        if n<L % no right neighbours for the rightmost column
            j2=j+y1(n)-y1(n+1);
            if 1<=j2  & j2<=ydiff(n+1)
                M(i1+j,i2+j2)=-1; %right neighbour
            end
            if 1<=j2-1 & j2-1<=ydiff(n+1)
                M(i1+j,i2+j2-1)=-1; %right-lower neighbour
            end
            if 1<=j2+1 & j2+1<=ydiff(n+1)
                M(i1+j,i2+j2+1)=-1;  %right-upper neighbour
            end
        end
        % compute element of data F and source G= Laplace(F)
        x=a+n*dx;
        y=(y1(n)+j)*dx;
        F(i1+j)= f(x,y); 
        G(i1+j)= ( f(x+dx,y)+f(x-dx,y)+f(x,y+dx)+f(x,y-dx)-4*f(x,y) )/dx^2;        
    end
    i1=i2;
end
M=M+M'+diag(8*ones(N,1),0); %fills in diagonal elements and left/lower neighbours in stiffness matrix

U=F+(M\G)*3*dx^2; %solves the main FEM equation U-F = (M/3)^(-1) (dx^2 G)

m=min(y1);
x=a+(1-5:L+5)*dx;
y=(m-5:max(y2)+5)*dx;

% re-order the N-vector U back as roughly an L x L matrix (two-variable function)  
Umatrix=(1.1* min(F)-0.1*max(F))* ones(L+10,max(y2)-min(y1)+1+10);  
i1=0; % start index of current column in the domain
for n=1:L
    for j=1:ydiff(n)
        Umatrix(n+5,y1(n)-m+j+5)=U(i1+j);
    end
    i1=i1+ydiff(n);
end

% plotting the solution U
figure(1)
clf
surf(x,y,Umatrix')
xlabel('x')
ylabel('y')
title('FEM solution')


%
% This is a Matlab script BEM.m which produces a 3d plot of the solution to the Dirichlet problem, which it calculates using the boundary element method with double layer potential and discretization with mesh size 1/N. It requires Matlab programs h.m for the boundary of the domain in polar coordinates and f.m for (an extension of) the boundary data
%
%


N= XXX ;  % mesh size = 1/ N
% Required functions: 
%   boundary curve in polar coordinates radius = h( angle )
%   extension of boundary Dirichlet data z=f(x,y)

% computes the mesh points p_i on the boundary of the domain
px= zeros(N,1);
py= zeros(N,1);
for i=1:N
  px(i)=h(2*pi*i/N)*cos(2*pi*i/N);
  py(i)=h(2*pi*i/N)*sin(2*pi*i/N);
end

d=zeros(N,1);  % length of boundary subset i
nx=zeros(N,1); % x-component of normal vector
ny=zeros(N,1);  % y-component of normal vector

% periodic definition of d, nx, ny
d(1)=sqrt((px(1)-px(N))^2+(py(1)-py(N))^2);
nx(1)= (py(1)-py(N))/d(1);
ny(1)= -(px(1)-px(N))/d(1);

for i=2:N
  d(i)=sqrt((px(i)-px(i-1))^2+(py(i)-py(i-1))^2);
  nx(i)= (py(i)-py(i-1))/d(i);
  ny(i)= -(px(i)-px(i-1))/d(i);
end

M=zeros(N,N);  % discretization of the kernel of the integral operator (multiplied by 2pi)
F=zeros(N,1);
for i=1:N
    M(i,i)=pi;  % diagonal elements
    F(i)= f(px(i),py(i));
    for j=1:i-1  % subdiagonal elements
        M(i,j)=d(j)*( (px(j)-px(i))*nx(j)+(py(j)-py(i))*ny(j) )/( (px(j)-px(i))^2+(py(j)-py(i))^2 );
    end
    for j=i+1:N  % superdiagonal elements
        M(i,j)=d(j)*( (px(j)-px(i))*nx(j)+(py(j)-py(i))*ny(j) )/( (px(j)-px(i))^2+(py(j)-py(i))^2 );
    end
end

V=M\F; %solves the main integral equation V= M^(-1) F

% below computes and plots the solution
dx=min(d); % mesh size to compute the solution at
x=(min(px)-5*dx):dx:max(px)+5*dx; % interior points of evaluations
y=(min(py)-5*dx):dx:max(py)+5*dx;
U=(1.1*min(F)-0.1*max(F))*ones(length(x), length(y));  % init of solution
for i=1:length(x)
    for j=1:length(y)
         if x(i)^2+y(j)^2<(h(angle(complex(x(i),y(j))))-dx/2)^2
            % checks if (x(i), y(j)) belongs to the domain, and is not too
            % near the boundary
            s=0;  % summation variable when computing the Riemann sum
            for k=1:N
                s=s+V(k)*d(k)*( (px(k)-x(i))*nx(k)+(py(k)-y(j))*ny(k) )/( (px(k)-x(i))^2+(py(k)-y(j))^2 );
            end
            U(i,j)=s;
        end
    end
end

% plot the solution U
figure(2)
clf
surf(x,y,U')
xlabel('x')
ylabel('y')
title('BEM solution')


%
% The required Matlab programs
%
%

function [ z ] = f( x,y )
  z= XXX ;
end

function [ y ] = h1( x )
  y= XXX ;
end

function [ y ] = h2( x )
  y= XXX ;
end

function [ r ] = h( v )
  r= XXX ;  
end