clear; clc; close all;
p = 70;
%p=-70;
lambda = @(x) 1; %constant diffusion
f = @(x) 0; %homogeneous RHS

N = 30;     % Initial mesh (uniform)
x = linspace(0,1,N)';

nIter = 8;

for it = 1:nIter
    [K,F] = assemble_1D(x,lambda,f); %assemble 
    for i = 1:length(x)-1 %convection term added manually
        h = x(i+1) - x(i);
        Ke_conv = p/2 * [-1 1; -1 1]; 
        K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + Ke_conv;
    end
    K(1,:) = 0; K(1,1) = 1; F(1) = 0; % Dirichlet BC at x=0
    K(end,:) = 0; K(end,end) = 1; F(end) = 1; % Dirichlet BC at x=1

    u = K\F;

    if it < nIter
        Nn = length(x); %initialize indicator
        eta = zeros(Nn,1);
        for j = 2:Nn-1 %interior nodes only
            ul = (u(j)-u(j-1))/(x(j)-x(j-1)); %left derivative
            ur = (u(j+1)-u(j))/(x(j+1)-x(j)); %right derivative
            %Jump in convection-diffusion flux
            Jl = -ul + p*u(j);
            Jr = -ur + p*u(j);
            % GH
            % JL = (u(j)-u(j-1))/(x(j)-x(j-1))*(x(j)-x(j-1));
            % JR = (u(j+1)-u(j))/(x(j+1)-x(j))*(x(j+1)-x(j));
            % HG
            eta(j) = abs(Jr - Jl);
            %GH eta(j) = abs(Jr - Jl)*sqrt((x(j)-x(j-1))*(x(j+1)-x(j)));
            %GH eta(j) = abs(Jr - Jl)*sqrt(max((x(j)-x(j-1)),(x(j+1)-x(j))));
        end
        x = r_adapt(x, eta);
        %
    end
end

figure
plot(x,u,'-o','LineWidth',1.5)
xlabel('x'), ylabel('u')
title(['Ex6C, r-adaptivity, Péclet p = ',num2str(p)])
grid on
% GH
hold on
fplot(@(x) (exp(p*x)-1)/(exp(p)-1),[0,1])
% Finalize the plot with the exact solution for comparison
legend('Numerical Solution', 'Exact Solution', 'Location', 'Best');