% Sheet 6 / Ex A
clf, clear, close all

% you need to play around with the parameters (enough starting nodes,
% refinements etc)

%p_vec = [5,10,20,100];
p = 100; % choose your p
n_nodes = 10;    % number of nodes on the coarse mesh (even - x=0 not included, odd - x=0 included)
nodes = linspace(-1,1,n_nodes);
f = @(x) 2*p^3.*x./(p^2.*x.^2+1).^2;
u_exact = @(x) atan(p*x);
lambda = @(x) 1+0*x;

% h-adaptivity
it_h = 4;   % number of refinements with h-adaptivity
for it = 1:it_h
    [K,f_vec] = assembling(nodes,lambda,f);
    % adaption for boundary (dirichlet left)
    f_vec(1) = f_vec(1) + K(1,1)*1e6*(-atan(p));
    K(1,1) = K(1,1)*(1+1e6);
    % adaption for boundary (neumann right)
    f_vec(end) = f_vec(end) + p/(p^2+1);
    
    % solving the system
    u = K\f_vec;
    
    % adapting the mesh
    if it < it_h
        nodes = adapt_h(nodes,u,lambda);
    end
end

figure(1)
hold on
plot(nodes,u,'-o')
fplot(u_exact,[-1,1])
title(['h-adapt: u_h, u for p=' num2str(p) ' after ' num2str(it_h) ' refinements, nstart=' num2str(n_nodes) ', nend=' num2str(length(nodes))]);
xlabel('x values');
ylabel('u');
legend('u_h', 'u');
grid on;
hold off
file_name_1 = ['ex_6_A_h-adap_p' num2str(p) '_ref' num2str(it_h) '_n_nodes' num2str(n_nodes) '.jpg'];
saveas(figure(1), file_name_1);


% r-adaptivity
n_nodes_r = 20;    % number of nodes on the coarse mesh (even - x=0 not included, odd - x=0 included)
nodes_r = linspace(-1,1,n_nodes_r);
it_r = 3;   % number of refinements with r-adaptivity
for it = 1:it_r
    [K,f_vec] = assembling(nodes_r,lambda,f);
    % adaption for boundary (dirichlet left)
    f_vec(1) = f_vec(1) + K(1,1)*1e6*(-atan(p));
    K(1,1) = K(1,1)*(1+1e6);
    % adaption for boundary (neumann right)
    f_vec(end) = f_vec(end) + p/(p^2+1);
    
    % solving the system
    u_r = K\f_vec;
    
    % adapting the mesh
    if it < it_r
        nodes_r = adapt_r(nodes_r,u_r,lambda);
    end
end

figure(2)
hold on
plot(nodes_r,u_r,'-o')
fplot(u_exact,[-1,1])
title(['r-adapt: u_h, u for p=' num2str(p) ' after ' num2str(it_r) ' refinements and n=' num2str(n_nodes_r)]);
xlabel('x values');
ylabel('u');
legend('u_h', 'u');
grid on;
hold off
file_name_2 = ['ex_6_A_r-adap_p' num2str(p) '_ref' num2str(it_r) '_n_nodes' num2str(n_nodes_r) '.jpg'];
saveas(figure(2), file_name_2);