function sec1_4_burgersFDM(jmax)
%
% example driver for Burgers' equation (with exact solution)
% in section 1.4. it calls movfdm() and interpolation_error().
%

%
% Copyright (C) 2010 Weizhang Huang and Robert D. Russell
% all rights reserved.
%
% This program is provided "as is", without warranty of any kind.
% Permission is hereby granted, free of charge, to use this program
% for personal, research, and education purposes. Distribution or use
% of this program for any commercial purpose is permissible
% only by direct arrangement with the copyright owner.
%

global epsilon;

  epsilon = 1.0e-4;

cpu0=clock;

   % job = 1 for solution and error at the final time
   % job = 2 for mesh trajectories

   job = 1;

   %jmax = 61;
   npde = 1;
   nn = npde*jmax;

   x=zeros(jmax,1);
   u=zeros(npde,jmax);
   
% call the moving mesh function
% monitor = 0 (fixed mesh), 1, 2, 3, 4, or 5
% mmpde = 4, 5, 6, 7 (mmpde = 7 ==> modified MMPDE5)
% alpha_type = 1, 2, or 3: (integral def, integral def with flooring, or alpha = constant)

   monitor = 3;
   mmpde = 7;
   alpha_type = 2;
   alpha = 1.0;
   reltol=1e-6; abstol=1e-4;
      
% generate initial adaptive mesh

   x=linspace(0, 1, jmax)';
   u(1,:)=zeros(1,jmax);
   xt=zeros(jmax,1);
   ut=zeros(npde,jmax);
   tspan = [0 0.2 0.4 0.6 0.8 1.0];
   [t,X,U]=movfdm(npde,jmax,tspan,x,xt,u,ut,@PDE_F2,@PDE_G2,@BC_L2,@BC_R2,...
                  [],monitor,[],[],[],mmpde,alpha_type,alpha); 

   fprintf(' --- End of generating initial mesh ---');

% define initial solution

   N=size(t,1);
   x=X(:,N);
   for j=1:jmax, u(1,j)=uexact(0,x(j)); end
   xt=zeros(jmax,1);
   ut=zeros(npde,jmax);

% call the moving mesh function
   
   if job==1 % for solution and error at the final time
   
      tspan = [0 0.2 0.4 0.6 0.8 1.0];
      [t,X,U]=movfdm(npde,jmax,tspan,x,xt,u,ut,@PDE_F,@PDE_G,@BC_L,@BC_R,...
                   [],monitor,reltol,abstol,[],mmpde,alpha_type,alpha);
      fprintf('\n cpu time used = %e \n', etime(clock,cpu0));
      % output the solution
      N=size(t,1);
	  figure
	  labels={};
      marks='+osdv^<>';
      lines1='--:-';
      lines2='-  .';
      colors='mkbgrc';
      for n=1:N
         u=U(:,:,n);
         x=X(:,n);
         hold on
         labels{n} = ['t = ' num2str(t(n)) ];
         style=['-' marks(1+rem(n,8))];
         plot(x,u(1,:)',style);%,'LineWidth',2);
      end
      hold off
      legend(labels{:});
      xlabel('x');
      ylabel('u');
      axis([0 1 -0.2 1.2]);
	  box on;
	  
	  % compute the error at t = T
	  errmax = interpolation_error(@uexact,jmax,npde,t(N),x,u,inf);
	  errL2 = interpolation_error(@uexact,jmax,npde,t(N),x,u);
	  errH1 = interpolation_error(@uexact,jmax,npde,t(N),x,u,2,1,1,@uexactDx);
	  fprintf(' (jmax,errmax,errL2,errH1) = %4d %e %e %e\n', jmax, errmax, errL2, errH1);
	  
   else % for trajectories
   
      tspan = [0 1.0];
      [t,X,U]=movfdm(npde,jmax,tspan,x,xt,u,ut,@PDE_F,@PDE_G,@BC_L,@BC_R,...
                   [],monitor,reltol,abstol,[],mmpde,alpha_type,alpha);
	  fprintf('\n cpu time used = %e \n', etime(clock,cpu0));
      plot(X',t);
	  axis([0 1 0 1]);
      xlabel('x')
      ylabel('t')
   end

fprintf('cpu time used = %e \n', etime(clock,cpu0));

% -----------------------------------------------------

function f=PDE_F(t,x,u,ux,ut)
  f(1,:) = ut(1,:);
  
function f=PDE_G(t,x,u,ux,ut)
  global epsilon;
  f(1,:) = epsilon*ux(1,:)-u(1,:).*u(1,:)*0.5;
  
function f=BC_L(t,x,u,ux,ut)
%  f(1) = u(1)-uexact(t,x);
  f(1) = ut(1)-uexactDt(t,x);
  
function f=BC_R(t,x,u,ux,ut)
%  f(1) = u(1)-uexact(t,x);
  f(1) = ut(1)-uexactDt(t,x);
  
% -----------------------------------------------------

function f=PDE_F2(t,x,u,ux,ut)
  N = size(x,2);
  for j=1:N
    f(1,j) = ut(1,j)-uexact(0,x(j));
  end
  
function f=PDE_G2(t,x,u,ux,ut)
  N = size(x,2);
  f(1,:) = zeros(1,N);
  
function f=BC_L2(t,x,u,ux,ut)
  f(1) = ut(1)-uexact(0,x);;
  
function f=BC_R2(t,x,u,ux,ut)
  f(1) = ut(1)-uexact(0,x);
  
% -----------------------------------------------------

function u=uexact(t,x)
   global epsilon;
   a1 = (-x+0.5-4.95*t)/(20.0*epsilon);
   a2 = (-x+0.5-0.75*t)/(4.0*epsilon);
   a3 = (-x+0.375)/(2.0*epsilon);
   a = max([a1 a2 a3]);
   r1 = exp((a1-a));
   r2 = exp((a2-a));
   r3 = exp((a3-a));
   u = (0.1*r1+0.5*r2+r3)/(r1+r2+r3);
   
% -----------------------------------------------------

function du=uexactDt(t,x)
   global epsilon;
      a1 = (-x+0.5-4.95*t)/(20.0*epsilon);
      a2 = (-x+0.5-0.75*t)/(4.0*epsilon);
      a3 = (-x+0.375)/(2.0*epsilon);
      a = max([a1, a2, a3]);
      r1 = exp(a1-a);
      r2 = exp(a2-a);
      r3 = exp(a3-a);
      r1t = -r1*4.95/(20*epsilon);
      r2t = -r2*0.75/(4*epsilon);
      r3t = 0.0;
      du = ((r1+r2+r3)*(0.1*r1t+0.5*r2t+r3t)-(r1t+r2t+r3t)*(0.1*r1+0.5*r2+r3)) ...
           /(r1+r2+r3)^2;
      
% -----------------------------------------------------

function du=uexactDx(t,x)
   global epsilon;
      a1 = (-x+0.5-4.95*t)/(20.0*epsilon);
      a2 = (-x+0.5-0.75*t)/(4.0*epsilon);
      a3 = (-x+0.375)/(2.0*epsilon);
      a = max([a1, a2, a3]);
      r1 = exp(a1-a);
      r2 = exp(a2-a);
      r3 = exp(a3-a);
      r1x = -r1/(20*epsilon);
      r2x = -r2/(4*epsilon);
      r3x = -r3/(2*epsilon);
      du = ((r1+r2+r3)*(0.1*r1x+0.5*r2x+r3x)-(r1x+r2x+r3x)*(0.1*r1+0.5*r2+r3)) ...
           /(r1+r2+r3)^2;

