% Solve a singularly perturbed reaction diffusion model problem in 1D
% An adaptive program with the residual error estimator of Verfuerth.
%
%	- eps*u'' + u = 0	in Omega = (0,1)
%	u(0)=1	  u(1)=0
%
% Analytical solution: Set eps2 := sqrt(eps)
% u = ( exp(-x/eps2) - exp((x-2)/eps2) ) / (1-exp(-2/eps2))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n0 = 10;	% number of initial intervals (equidistant mesh)
n = n0;		% present number of intervals
x = 0:1/n:1;	% the initial mesh

gamma = 0.1;	% Refinement occurs if eta_T > gamma*eta_max
		% eta_T := error estimator, eta_max := max_T {eta_T}

eps=input('Please enter epsilon:          ');
if isempty(eps) == 1 
	eps = 1e-6;
	fprintf('Use default value for epsilon: %4.1e\n',eps);
end
eps2 = sqrt(eps);
eps32 = eps2^3;

fprintf('Number of initial intervals:   %2d\n',n0); 
fprintf('Refinement parameter gamma:   %5.2f\n\n',gamma);

fprintf('                      Upper Bd  Lower Bd\n');
fprintf('Elements: |||Err|||   Err/Est   Est/Err     C-Err   Conv-rate\n');
fprintf('------------------------------------------------------------\n');

fid = fopen('result-Enorm','w');
fprintf(fid,'# Energienorm, Epsilon =   %4.1e\n\n',eps);
fclose(fid);

fid = fopen('result-Cnorm','w');
fprintf(fid,'# Maximum Norm, Epsilon =   %4.1e\n\n',eps);
fclose(fid);

while 1		% Adaptive loop


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the Finite Element Solution.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h = x(2:n+1) - x(1:n);	% mesh sizes

D = [1 , eps*( 1./h(1:n-1) + 1./h(2:n) ) , 1]' + ( [h 0] + [0 h])'/3;
NDu = [-eps*1./h + h/6 , 0]';
NDo = [0 , -eps*1./h + h/6 ]';

A = spdiags([NDu D NDo] , -1:1 , n+1,n+1);	% stiffness matrix
A(1,1) = 1;
A(1,2) = 0;
A(n+1,n) = 0;
A(n+1,n+1) = 1;

b = [1; zeros(n,1) ];	% right hand side
uh = A\b;		% FEM solution

% Exact solution: Numerically stable for small eps
u = ( exp(-x'/eps2) - exp((x'-2)/eps2) ) / (1-exp(-2/eps2));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the error estimator.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The local error estimator eta_T is for an interval T with two
% nodes E_i.
% The face residual r_E is the canonical choice, i.e. eps*Gradient jump.
%
% eta_T^2 := alpha_T^2 * || r_T ||_T^2 +
%	     + eps^{-1/2}*alpha_T * Summe_E_i | r_E |^2
% with alpha_T := min(1,h_T/sqrt(eps))

al = min( 1 , h / sqrt(eps))';

% Element residual:
norm_rT2 = h'/3 .* ( uh(1:n).^2 + uh(2:n+1).^2 + uh(1:n).*uh(2:n+1) );

% "Face" residual = gradient jump * eps. Zero at x=0 and x=1.
rE = eps*[ 0 ;( uh(3:n+1) -uh(2:n)) ./h(2:n)' - ( uh(2:n) -uh(1:n-1)) ./h(1:n-1)' ; 0] ;

% Error estimator (squared)
Est = al.^2 .* norm_rT2 + eps^(-1/2) * al .* (rE(1:n).^2 + rE(2:n+1).^2);

Est_glob = sqrt(sum(Est));	% global error estimator
max_err = max(abs(uh - u));	% L_{infinity} error
Est_max = max(Est);		% Maximum local error estimator


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the (analytical) error in the energy norm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1 = x(2:n+1)';	% Auxiliary (column) vector x_{i+1}
x0 = x(1:n)';

% Error fct = [ e^{-x/eps2} - e^((x-2)/eps2) ] / (1-e^{-2/eps2}) - c*x-d
% ||| error |||_T gives a formula of about one page length (obtained via
% MAPLE).
% For small eps this formula in unstable (gives NaN). Then one has
% approximately u ~ e^{-x/eps2}, and the error is (approximately and
% very accurately) as below.

c = ( uh(2:n+1) - uh(1:n) ) ./ h';
d = uh(1:n) - c.*x0;

if eps >= 0.001 	% analytical formula by maple

H0 = -1/6*(-12.*c.*x0.^2.*d.*exp(2/eps2)-6.*c.*x1.^2.*d.*exp(4/eps2)+6.*c.*x0.^2.*d.*exp(4/eps2)+4.*c.^2.*x1.^3.*exp(2/eps2)+2.*c.^2.*x0.^3.*exp(4/eps2)-4.*c.^2.*x0.^3.*exp(2/eps2)+12*c*eps.*exp(-(x0-4)/eps2)-12.*c*eps.*exp(-(x0-2)/eps2)+6.*d.^2.*x0.*exp(4/eps2)-12.*d.^2.*x0.*exp(2/eps2)-12.*c*eps.*exp((x0+2)/eps2)+12.*c*eps2.*x0.*exp((x0+2)/eps2)+2.*c.^2.*x0.^3+6.*d.^2.*x0+3*eps2.*exp(-2*(x1-2)/eps2)+12.*x1.*exp(2/eps2)+3*eps2.*exp(2.*x0/eps2)-12.*x0.*exp(2/eps2)-3*eps2.*exp(2.*x1/eps2)+12.*d*eps2.*exp((x0+2)/eps2)-2.*c.^2.*x1.^3.*exp(4/eps2)-12.*d*eps2.*exp(-(x0-2)/eps2)+12.*c*eps.*exp(x0/eps2)+12.*d*eps2.*exp(-(x0-4)/eps2)-3*eps2.*exp(-2*(x0-2)/eps2)+6.*c.*x0.^2.*d-6.*c.*x1.^2.*d+12.*c.*x1.^2.*d.*exp(2/eps2)-2.*c.^2.*x1.^3-6.*d.^2.*x1-12.*c*eps2.*x1.*exp(-(x1-4)/eps2)+12.*c*eps2.*x1.*exp(-(x1-2)/eps2)+12.*c*eps2.*x0.*exp(-(x0-4)/eps2)-12.*c*eps2.*x0.*exp(-(x0-2)/eps2)-12.*c*eps2.*x0.*exp(x0/eps2)+12.*c*eps2.*x1.*exp(x1/eps2)-12.*c*eps2.*x1.*exp((x1+2)/eps2)-12.*d*eps2.*exp(x0/eps2)+12.*c*eps.*exp(-(x1-2)/eps2)+12.*c*eps.*exp((x1+2)/eps2)-12.*d*eps2.*exp((x1+2)/eps2)-12.*d*eps2.*exp(-(x1-4)/eps2)-6.*d.^2.*x1.*exp(4/eps2)+12.*d.^2.*x1.*exp(2/eps2)+12.*d*eps2.*exp(-(x1-2)/eps2)+12.*d*eps2.*exp(x1/eps2)-12.*c*eps.*exp(-(x1-4)/eps2)-12.*c*eps.*exp(x1/eps2))/(exp(2/eps2)-1).^2;


H1 = 1/2*(-eps.*exp(-2*(x1-2)/eps2)+4.*x1.*exp(2/eps2)*eps2-4*c*eps32.*exp(-(x1-4)/eps2)+4*c*eps32.*exp(-(x1-2)/eps2)+exp(2.*x1/eps2)*eps+4*c*eps32.*exp((x1+2)/eps2)-4*c*eps32.*exp(x1/eps2)+2*c.^2.*x1*eps32.*exp(4/eps2)-4*c.^2.*x1*eps32.*exp(2/eps2)+2*c.^2.*x1*eps32+eps.*exp(-2*(x0-2)/eps2)-4.*x0.*exp(2/eps2)*eps2+4*c*eps32.*exp(-(x0-4)/eps2)-4*c*eps32.*exp(-(x0-2)/eps2)-exp(2.*x0/eps2)*eps-4*c*eps32.*exp((x0+2)/eps2)+4*c*eps32.*exp(x0/eps2)-2*c.^2.*x0*eps32.*exp(4/eps2)+4*c.^2.*x0*eps32.*exp(2/eps2)-2*c.^2.*x0*eps32)/eps32/(exp(2/eps2)-1).^2;


else	% approximate formula for small epsilon

    % First  Int (err')^2
    H1 = ( exp(-x1/eps2)-exp(-x0/eps2) ) .* ...
	( ( exp(-x1/eps2)+exp(-x0/eps2) )/(-2)/eps2 - 2*c) ...
	+ c.^2.*h';

    % Now Int (err)^2
    H0_a = - eps2/2 * ( exp(-2*x1/eps2)-exp(-2*x0/eps2) );
    H0_b = 2*eps2* exp(-x1/eps2) .* (c.*x1 + c*eps2 + d) - ...
      (2*eps2* exp(-x0/eps2) .* (c.*x0 + c*eps2 + d) );
    H0_c = c.^2/3.*(x1.^3-x0.^3) + c.*d.*(x1.^2-x0.^2) + d.^2 .*h';
    H0 = H0_a + H0_b + H0_c;

end	% end if

% The element error (squared) in the energy norm. Verified with MAPLE.
Err = eps*H1 + H0;	

Err_glob = sqrt(sum(Err));	% global error
Err_max = max(Err);		% Maximum local error, squared

% The ratio of the lower error bound. Has to be bounded from above.
low_bd = Est ./ ( Err + [0 ;Err(1:n-1)] + [Err(2:n) ;0] );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The convergence rate (energy norm) between two successive steps.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
err_new = Err_glob;
n_new = n;

if n > n0
	konv = -log(err_new/err_old) / log(n_new/n_old);
else
	konv = 0;
end;

n_old=n_new;
err_old=err_new;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The result and some nice plots:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,1); plot(x,[uh  u]);
title('u and u_h')

subplot(2,2,2); plot(x,uh - u);
title('u - u_h');

subplot(2,2,3); plot(x,log10([h 0.1001]));
axis([0 1 log10(h(1))-1 0])
title('logarithmic mesh density');

subplot(2,2,4); plot(x,[sqrt(low_bd); sqrt(low_bd(n))]);
title('local Ratio Lower Bound: Est/Err');

fprintf('  %5d    %7.2e  %7.4f   %7.4f    %8.3e   %6.2f\n',...
n, Err_glob,Err_glob/Est_glob,max(sqrt(low_bd)),max_err,konv);

% Write to some file "result...".
fid = fopen('result-Enorm','a');
fprintf(fid,'  %5d    %10.4e \n',n, Err_glob);
fclose(fid);

fid = fopen('result-Cnorm','a');
fprintf(fid,'  %5d    %10.4e \n',n, max_err);
fclose(fid);

pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Refinement:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Loop over all elements:
% When Eta_T > gamma * Eta_max then insert new node.

m = 1;
clear y
for i=1:n
	y(m) = x(i);
	if Est(i) > gamma^2*Est_max	% mesh control via Estimator
%	if Err(i) > gamma^2*Err_max	% mesh control via true error
		y(m+1) = (x(i)+x(i+1))/2;
		m = m+2;
	else
		m = m+1;
	end;
end;
y(m) = 1;

x = y;
n = m-1;

end	% End while