
% parameter of ODE 
b = @(x) ones(size(x,1),1);
%b = @(x) x;
%b = @(x) -(x-0.5);
%b = @(x) cos(2*pi*x);
c = @(x) zeros(size(x,1),1);
f = @(x) ones(size(x,1),1);

epsilon = 0.001;

% exact solution for b=1, c=0, f=1
uexact = @(x) x-(exp(-(1-x)/epsilon)-exp(-1/epsilon))/(1-exp(-1/epsilon));


% boundary conditions 
bc = [1,0,0; 1,0,0]; 

% parameters for discretisation
hunif=0.5;
Nref = 15;
error = zeros(Nref,1);
ndof4lvl = zeros(Nref,1);
mesh_gen = @(x) mesh_uniform(x);

% Bakhvalov
% A=1;
% q=0.5;
% mesh_gen = @(x) bakhvalov(x,epsilon,A,q);

% Shishkin 
% beta = 1;
% mesh_gen = @(x) Shishkin(x,epsilon,beta);

%compute reference solution 
h_refer = min(epsilon/4,2^(-Nref-2));
mesh_refer_unif = (0:h_refer:1)';
mesh_refer = mesh_gen(mesh_refer_unif);
u_refer = computeFDsolution(epsilon,b,c,f,bc,mesh_refer,@(q) ones(size(q,1),1));

for j = 1:Nref
    % parameter for discretization
    hunif = hunif/2;
    mesh_unif = (0:hunif:1)';
    mesh = mesh_gen(mesh_unif);
    
    % non uniform mesh
    %mesh = [(0:0.1:0.9)';(0.91:0.01:1)'];
    
    sigma_fct = @(q) ones(size(q,1),1); % no artificial diffusion
    %sigma_fct = @(q) 1 + q;  % simple upwind scheme
    %sigma_fct = @(q) max(1,q); 
    %sigma_fct = @(q) sqrt(1+q.^2);
    %sigma_fct = @(q) 1+q.^2./(1+q); % Samarskii's scheme
    %sigma_fct = @(q) q.*coth(q); % Il'in-Allen-Southwell scheme
    % SD-FEM for P1, if c=0 (up to quadrature) for the choice from the
    %     lecture, for uniform meshes
    %    modification of rhs has to be done!
    % sigma_fct = @(q) 1+0.2*q   % if epsilon<h
    % sigma_fct = @(q) 1+0.2*hunif/epsilon   % if h<epsilon

    u = computeFDsolution(epsilon,b,c,f,bc,mesh,sigma_fct);
    
    
    % compute error with respect to reference solution
    uref_fct = @(x) evaluate_uref(mesh_refer,u_refer,x); 
    % compute l^infty error 
    error(j) = max(abs(u-uref_fct(mesh)));
    %error(j) = max(abs(u-uexact(mesh)));

    ndof4lvl(j) = size(mesh,1)-2; 

    % figure; 
    % plot(mesh,u,'-x');
    % hold all;
    % plot(mesh,zeros(size(mesh,1),1),'x');

end

figure;
loglog(ndof4lvl,error);


% % plot first derivative
% du = (u(2:end)-u(1:(end-1)))./h;
% figure;
% plot((mesh(1:(end-1))+mesh(2:end))/2,du);
% 
% % plot second derivative
% d2u = -A*u;
% figure;
% plot(mesh(2:(end-1)),d2u);



function val = evaluate_uref(mesh,v,x)
% computes v(x) where v is a piecewise linear function on mesh
% only works for sorted x (implementation is then faster)
nnodes = size(mesh,1);

xj = zeros(size(x,1),1);
jx = zeros(size(x,1),1);

j=1;
for k = 1:size(x)
    diff = x(k)-mesh; 
    
    for ell=j:nnodes % seek largest j=j(ell) such that mesh(j)<=x(ell)        
        if diff(ell) <=0
            break;
        end
        j = j+1;
    end
    jx(k) = j;
    xj(k) = mesh(j);
end

% if jx(k)=nnodes, then x(k)=1; 
% in this case, it doesn't matter, if xj=1 or xjp1=1
ind1 = find(jx==nnodes);
jx(ind1) = jx(ind1)-1;
xj(ind1) = mesh(jx(ind1));

xjp1 = mesh(jx+1);
val = v(jx) + (x-xj).*(v(jx+1)-v(jx))./(xjp1-xj);

end

function val = mesh_uniform(mesh)
% uniform mesh 
val = mesh;
end

function mesh_new = bakhvalov(mesh,epsilon,A,q)
% generates Bakhvalov mesh 

if size(mesh,1)>500
    mesh;
end

%compute tau 
tau = q-A*epsilon;
for k=1:100
    tau = q-(A*epsilon*(1-tau))/(1+A*epsilon*log(1-tau/q));
end

ind1 = find(mesh>1-tau);

psi1minustau = 1+A*epsilon*log(1-(tau/q));
dpsi1minustau = A*epsilon/(q-tau);
mesh_new = psi1minustau+(mesh-1+tau)*dpsi1minustau;
mesh_new(ind1) = 1+A*epsilon*log(1-(1-mesh(ind1))/q);
mesh_new(1) = 0; % corrects the approximation error of tau
mesh_new(end) = 1; % to avoid round off errors at 1
end

function mesh_new = Shishkin(mesh,epsilon,beta)
% generates Shishkin mesh 

sigma = (2/beta)*epsilon*log(size(mesh,1)-1);

if sigma>0.5
    mesh_new = mesh;
else 
    mesh_new = 2*(1-sigma)*mesh;
    ind1 = find(mesh>1/2);
    % alternative, if mesh is uniform:
    % ind1 = ceil((size(mesh,1)+1)/2):size(mesh,1);
    mesh_new(ind1) = 2*sigma*mesh(ind1)+1-2*sigma;
    mesh_new(end) = 1; % to avoid round off errors at 1
end

end

