function [t,x,u]=movfem(npde,jmax,tspan,x0,xt0,u0,ut0,PDE_F,PDE_G,PDE_Fa,PDE_Fb,...
                        BC_L,BC_R,TauFun,monitor,reltol,abstol,MonitorFun,mmpde,alpha_type,alpha)
%
% [t,x,u]=movfem(npde,jmax,tspan,x0,xt0,u0,ut0,PDE_F,PDE_G,PDE_Fa,PDE_Fb,...
%                BC_L,BC_R,TauFun,monitor,reltol,abstol,MonitorFun,mmpde,alpha_type,alpha)
%
% This function integrates a system of (npde) one-dimensional PDEs
% in the weak form
%
%   int_a^b ( F(t,x,u,u_x,u_t)*Phi + G(t,x,u,u_x,u_t)*Phi_x ) dx
%      + Fa(t,x,u,u_t)*Phi|_(x=a) + Fb(t,x,u,u_t)*Phi|_(x=b)
%                       = 0    for all possible test functions Phi
%                      0 = B_L(t,x,u,u_t)         at x = a
%                      0 = B_R(t,x,u,u_t)         at x = b
%
% using an MMPDE moving mesh strategy. The modified MMPDE5 is used as
% the default strategy.
%  
% The physical PDEs are discretized using the linear finite element method
% and mesh PDEs discretized using central finite differences.
% The resulting ODE system is integrated using the Matlab built-in solver ode15i.
%
% Note that the first 13 arguments must be specified by the user. The other 8
% arguments (monitor,reltol,abstol,mmpde,alpha_type,alpha) are optional.
%
% output variables:
%
% t:        column vector of time instants, t(1:N)
% x:        meshes at all time levels, x(j,n), j=1:jmax, n=1:N
% u:        solutions at all time levels, u(i,j,n), i=1:npde, j=1:jmax, n=1:N
%
% input variables:
%
% npde:     the number of physical PDEs contained in the system to be integrated.
% jmax:     the number of mesh points.
% tspan:    an array specifying the time interval for integration.
% x0:       column vector of size jmax, initial mesh.
% xt0:      column vector of size jmax, approximate values for initial mesh speed
%           which are required by the implicit DAE solver ode15i.
% u0:       array of size npde by jmax, intial values for the physical solution.
%           u0(i,j): i=1:npde, j=1:jmax
% vt0:      array of size npde by jmax, approximate values for the initial time
%           derivatives of the physical solution which are required by the implicit
%           DAE solver ode15i.
% PDE_F,PDE_G,PDE_Fa,PDE_Fb,BC_L,BC_R: functions used to define the system of
%           physical PDEs:
%
%            PDE_F := F(t,x,u,u_x,u_t)
%            PDE_G := G(t,x,u,u_x,u_t)
%           PDE_Fa := Fa(t,x,u,u_t)
%           PDE_Fb := Fb(t,x,u,u_t)
%             BC_L := B_L(t,x,u,u_t)
%             BC_R := B_R(t,x,u,u_t)
%
%    function f=PDE_F(t,x,u,ux,ut)
%    function g=PDE_G(t,x,u,ux,ut)
%
%           input variables:
%           t is scalar, x=x(1:m) is a row vector, and u=u(1:npde,1:m), ux=ux(1:npde,1:m),
%           and ut=ut(1:npde,1:m) are matrices.
%
%           output variable:
%           f=f(1:npde,1:m) and g=g(1:npde,1:m) should be matrices.
%
%    function res=PDE_Fa(t,x,u,ut)
%    function res=PDE_Fb(t,x,u,ut)
%    function res=BC_L(t,x,u,ut)
%    function res=BC_R(t,x,u,ut)
%
%           input variables:
%           t and x are scalar and u=u(1:npde), ux=ux(1:npde), and ut=ut(1:npde) are
%           column vectors.
%
%           output variable:
%           res=res(1:npde) should be a column vector.
%
%    note: BC_L and BC_R are used to define the dirichlet bcs only. when, for exmaple,
%           the l-th component has a neumann or mixed bc at x=a, then set BC_L(l) = NaN.
%
% TauFun:   the function defines the parameter tau. it should have the form
%
%           function tau=TauFun(npde,jmax,t,x,u)
%
%           where (npde,jmax,t,u) are input variable and u contains the values of
%           the physical variables and the mesh at time t. x is a column vector
%           of size jmax and u is a matrix of size npde by jmax: u(i,j), i=1:npde,
%           j=1:jmax.
%
%           default value tau = 0.01 is invoked when TauFun is not defined
%           thru absence or putting [] for this argument.
%
% monitor:  an integer specifying the mesh density function.
%           0: for fixed mesh.
%           1: k=0, m=0 (piecewise constant approximation, error measured in L2).
%           2: k=1, m=0 (piecewise linear approximation, error measured in L2).
%           3: k=1, m=1 (piecewise linear approximation, error measured in
%                        H1 semi-norm).
%           4: arclength mesh density function, rho = sqrt(1+(u')^2).
%           5: curvature mesh density function, rho = (1+(u'')^2)^0.25.
%
%           default: monitor=3
%
% reltol, abstol: relative and absolute tolerances used in integration of
%           the ode system.
% MonitorFun: the user defined monitor function. it should have the form
%
%           function rho=MonitorFun(npde,jmax,t,x,u)
%
%           where (npde,jmax,t,x,u) are input variable and u contains the values of
%           the physical variables and the mesh at time t. x is a column vector
%           of size jmax and u is a matrix of size npde by jmax: u(i,j), i=1:npde,
%           j=1:jmax.
%
%           an internal monitor function (specified through (monitor)) is invoked
%           when MonitorFun is not defined thru absence or putting [] for
%           this argument.
%
%   note: when the user provides his own montor function, the function must return
%           a positive, column vector of size (jmax). in this case, the variable
%           (monitor) should be set to a value other than 0.
% mmpde:    (4, 5, 6, 7) specifies which MMPDE (out of MMPDE4, MMPDE5, MMPDE6,
%           modified MMPDE5) is used. these MMPDEs have the expressions as
%
%           MMPDE5: x_t = 1/(tau) (rho x_xi)_xi,
%           MMPDE4: (rho (x_t)_xi)_xi = -1/(tau) (rho x_xi)_xi,
%           MMPDE6: ((x_t)_xi)_xi = -1/(tau) (rho x_xi)_xi,
%  modified MMPDE5: x_t = 1/(rho*tau) (rho x_xi)_xi,  (default)
%
%           where rho is the monitor function (or more precisely, the mesh density
%           function) and tau, which can be defined as a function of t, mesh,
%           and the physical variables, is a parameter for adjusting the response
%           of the mesh to the changes in rho.
% alpha_type and alpha: these two variables are used to control the definition of
%           the adaptation intensity parameter alpha. alpha_type indicates
%           the type of the definition.
%           alpha_type = 1: integral, solution based definition. in this case,
%                           the value of variable (alpha) is not used.
%           alpha_type = 2: integral, solution based definition with flooring. that is,
%                           alpha := max(alpha, integral definition).
%           alpha_type = 3: alpha := alpha (constant).
%
%           default: alpha_type = 1, alpha = 1
%

%
% 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.
%

% functions contained or called in this file:
%
% 1. res=movfdm_res(t,v,vt,PDE_F,PDE_G,BC_L,BC_R,jmax,npde,TauFun,monitor)
% 2. [SPDY,SPDYP] = jpattern(npde,jmax,mmpde,alpha_type)
% 3. rho = mesh_density_fnct(npde,jmax,x,v,monitor,lsfitting,alpha_type,alpha) 
%                                             (called external function)
% 4. dv = grad_recovery(jmax,x,v)             (called external function)
% 5. rho1 = mesh_density_fnct_smoothing(jmax,rho)      (called external function)
% 6. tau = TauFun_default(npde,jmax,t,x,v)
%

   % set the default values for the parameters
   
   if nargin < 13
       error('The first 13 arguments must be defined.');
   elseif nargin < 14
      TauFun = @TauFun_default;
      monitor=[];
      reltol=[];
      abstol=[];
      MonitorFun=[];
	  mmpde=[];
	  alpha_type=[];
	  alpha=[];
   elseif nargin < 15
      monitor=[];
      reltol=[];
      abstol=[];
      MonitorFun=[];
	  mmpde=[];
	  alpha_type=[];
	  alpha=[];
   elseif nargin < 16
      reltol=[];
      abstol=[];
      MonitorFun=[];
	  mmpde=[];
	  alpha_type=[];
	  alpha=[];
   elseif nargin < 17
      abstol=[];
      MonitorFun=[];
	  mmpde=[];
	  alpha_type=[];
	  alpha=[];
   elseif nargin < 18
      MonitorFun=[];
	  mmpde=[];
	  alpha_type=[];
	  alpha=[];
   elseif nargin < 19
	  mmpde=[];
	  alpha_type=[];
	  alpha=[];
   elseif nargin < 20
	  alpha_type=[];
	  alpha=[];
   elseif nargin < 21
	  alpha=[];
   elseif nargin > 21
      error('Two many arguments are specified.');
   end    

   flag_TauFun = 0;
   flag_monitor = 0;
   flag_reltol = 0;
   flag_abstol = 0;
   flag_mmpde = 0;
   flag_alpha_type = 0;
   flag_alpha = 0;
   if isempty(TauFun)
      TauFun = @TauFun_default;
      flag_TauFun = 1;
   end
   if isempty(monitor) | (monitor<0) | (monitor>5)
      monitor=3;
      flag_monitor = 1;
   end
   if isempty(reltol) | (reltol<=0)
      reltol=1e-6;
      flag_reltol = 1;
   end
   if isempty(abstol) | (abstol<=0)
      abstol=1e-4;
      flag_abstol=1;
   end
   if isempty(mmpde) | (mmpde<4) | (mmpde>7)
      mmpde=7;
	  flag_mmpde = 1;
   end
   if isempty(alpha_type) | (alpha_type<1) | (alpha_type>3)
	  alpha_type=1;
	  flag_alpha_type = 1;
   end
   if isempty(alpha) | (alpha<=eps)
      alpha = 1.0;
	  flag_alpha = 1;
   end

% set options

   nn = npde*jmax;

   [SPDY,SPDYP]=jpattern(npde,jmax,mmpde,alpha_type);
   j=round(jmax/2);
   if j < 1, j = 1; end
   j = npde*(j-1)+1;
   opts=odeset('RelTol',reltol,'AbsTol',abstol,'OutputFcn',@odeprint,...
               'OutputSel',j,'JPattern',{SPDY,SPDYP},'Stats','on');
      
% use a helper to find better initial derivatives

   free_v=zeros((npde+1)*jmax,1);
   free_v(nn+1)=1;
   free_v(nn+npde)=1;
   v0 =[reshape(u0,npde*jmax,1)
        x0];
   vt0 =[reshape(ut0,npde*jmax,1)
        x0];
   
   [v0md,vt0]=decic(@movefem_res,tspan(1),v0,free_v,vt0,[],opts,...
                    PDE_F,PDE_G,PDE_Fa,PDE_Fb,BC_L,BC_R,...
                    jmax,npde,TauFun,monitor,MonitorFun,mmpde,alpha_type,alpha);
   dd1 = norm(v0md(1:nn)-v0(1:nn),inf);
   dd2 = norm(v0md((nn+1):end)-v0((nn+1):end),inf);
   
% call ode15i

   [t,v]=ode15i(@movefem_res,tspan,v0md,vt0,opts,...
                PDE_F,PDE_G,PDE_Fa,PDE_Fb,BC_L,BC_R,...
                jmax,npde,TauFun,monitor,MonitorFun,mmpde,alpha_type,alpha);

   N=size(t,1);
   x=zeros(jmax,N);
   u=zeros(npde,jmax,N);
   for n=1:N
      u(:,:,n)=reshape(v(n,1:nn),npde,jmax);
      x(:,n)=v(n,(nn+1):end)';
   end

% print messages

if flag_TauFun==1
   fprintf('\n *** The default value tau=0.01 is used.\n');
end
if flag_monitor==1
   fprintf('\n *** The default value monitor=3 is used.\n');
end
if flag_reltol==1
   fprintf('\n *** The default value reltol=1e-6 is used.\n');
end
if flag_abstol==1
   fprintf('\n *** The default value abstol=1e-4 is used.\n');
end
if flag_mmpde==1
   fprintf('\n *** The default mmpde=7 (modified MMPDE5) is used.\n');
end
if flag_alpha_type==1
   fprintf('\n *** The default definition (integral definition) for alpha is used.\n');
end
if flag_alpha==1
   fprintf('\n *** The default alpha = 1.0 is used.\n');
end
if (dd1+dd2)>eps
   fprintf('\n *** Warning: a modified initial solution is used.\n');
     fprintf('     max diff in solution = %e\n', dd1);
     fprintf('         max diff in mesh = %e\n', dd2);
end

% end of movefem

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

function res=movefem_res(t,v,vt,PDE_F,PDE_G,PDE_Fa,PDE_Fb,BC_L,BC_R, ...
					     jmax,npde,TauFun,monitor,MonitorFun,mmpde,alpha_type,alpha)
%
% this function computes the res function to be called by an ode
% solver. the res function results from the linear finite element
% discretizetion of the (npde) physical pdes, the central finite
% difference discretization of the mesh equation, and discretization
% of the corresponding boundary conditions.
%
% mmpde = 4, 5, 6, or 7 (modified MMPDE5)
%

   % declare arrays
   
   Res=zeros(npde,jmax);
   Resx=zeros(jmax,1);
   
   u=zeros(npde,1);
   ux=zeros(npde,1);
   ut=zeros(npde,1);
   f=zeros(npde,1);
   g=zeros(npde,1);

   XX=zeros(1,2*(jmax-1));
   U=zeros(npde,2*(jmax-1));
   UX=zeros(npde,2*(jmax-1));
   UT=zeros(npde,2*(jmax-1));
   F=zeros(npde,2*(jmax-1));
   G=zeros(npde,2*(jmax-1));   
   rho=zeros(jmax,1);   
   nn = npde*jmax;

   npts = 2;   
   s = zeros(npts,1);
   w = zeros(npts,1);   
   s(1) = 0.5*(1-1.0/sqrt(3));
   s(2) = 0.5*(1+1.0/sqrt(3));
   w(1) = 0.5;
   w(2) = 0.5;
   
   % extract physical variables and mesh from v and vt

   V = reshape(v(1:nn),npde,jmax);
   Vt = reshape(vt(1:nn),npde,jmax);
   x = v((nn+1):end);
   xt = vt((nn+1):end);
   
   % compute monitor function and tau
   
   if monitor~=0
      if isempty(MonitorFun) % use internal monitor functions
	     lsfitting = 'no'; % using finite difference approximation
	     rho = mesh_density_fnct(npde,jmax,x,V,monitor,lsfitting,alpha_type,alpha);
      else
         rho = feval(MonitorFun,npde,jmax,t,x,V);
      end      
      % smooth the monitor function   
      for j=1:4
         Resx = rho;
         rho = mesh_density_fnct_smoothing(jmax,Resx);
      end
      % compute tau   
      tau = feval(TauFun,npde,jmax,t,x,V);
      if isempty(tau) | (tau<=0)
         error('tau=TauFun(npde,jmax,t,x,v) should be positive.');
      end
   end
      
   % for physical pdes: interior points

   for j=2:jmax
      % perform numerical quedrature on [x_{j-1},x_j]
      dx = x(j)-x(j-1);
	  dxinv = 1.0/dx;
      dxt = xt(j)-xt(j-1);      
      for i=1:npts
	     k = 2*(j-2)+i;
         xx = x(j-1)+s(i)*dx;
         xxt = xt(j-1)+s(i)*dxt;
         phi1 = (x(j)-xx)*dxinv;
         phi2 = (xx-x(j-1))*dxinv;
         phi1dx = -dxinv;
         phi2dx = dxinv;
		 XX(k) = xx;
         U(:,k) = V(:,j-1)*phi1 + V(:,j)*phi2;
         UX(:,k) = V(:,j-1)*phi1dx + V(:,j)*phi2dx;
         UT(:,k) = Vt(:,j-1)*phi1 + Vt(:,j)*phi2 - UX(:,k)*xxt;
      end
   end
   F= feval(PDE_F,t,XX,U,UX,UT);
   G= feval(PDE_G,t,XX,U,UX,UT);
   for j=2:jmax
      % perform numerical quedrature on [x_{j-1},x_j]
      dx = x(j)-x(j-1);
	  dxinv = 1.0/dx;    
      for i=1:npts
	     k = 2*(j-2)+i;
         xx = x(j-1)+s(i)*dx;
         phi1 = (x(j)-xx)*dxinv;
         phi2 = (xx-x(j-1))*dxinv;
         phi1dx = -dxinv;
         phi2dx = dxinv;
         Res(:,j-1) = Res(:,j-1)+dx*w(i)*(phi1*F(:,k)+phi1dx*G(:,k));
         Res(:,j) = Res(:,j)+dx*w(i)*(phi2*F(:,k)+phi2dx*G(:,k));
      end
   end
   
   % for physical pdes: boundary points

   j=1;   
   xx = x(j);
   xxt = xt(j);
   u = V(:,j);
   ut = Vt(:,j);
   f = feval(PDE_Fa,t,xx,u,ut);
   Res(:,j) = Res(:,j) + f;   % for neumann bcs
   f = feval(BC_L,t,xx,u,ut); % for dirichlet bcs
   for l=1:npde
%      if ~isempty(f(l)), Res(l,j) = f(l); end
       if ~isnan(f(l)), Res(l,j) = f(l); end
   end
         
   j=jmax;
   xx = x(j);
   xxt = xt(j);
   u = V(:,j);
   ut = Vt(:,j);
   f = feval(PDE_Fb,t,xx,u,ut);
   Res(:,j) = Res(:,j) + f;   % for neumann bcs
   f = feval(BC_R,t,xx,u,ut); % for dirichlet bcs
   for l=1:npde
%      if ~isempty(f(l)), Res(l,j) = f(l); end
	  if ~isnan(f(l)), Res(l,j) = f(l); end
   end

   % for mesh equation
   
   if monitor==0 % for fixed meshes

       Resx = xt;
    
   else  % for moving meshes
      
      % for mesh pde: interior points (mmpde5 with balancing factor)
      dxi = jmax-1;
	  if mmpde==4
         for j=2:(jmax-1)
            Resx(j) = (tau*(xt(j+1)-xt(j))+(x(j+1)-x(j)))*(rho(j+1)+rho(j)) ...
              -(tau*(xt(j)-xt(j-1))+(x(j)-x(j-1)))*(rho(j)+rho(j-1));
         end
      elseif mmpde==6
         for j=2:(jmax-1)
            Resx(j) = (tau*(xt(j+1)-xt(j))+(x(j+1)-x(j))*(rho(j+1)+rho(j))) ...
              -(tau*(xt(j)-xt(j-1))+(x(j)-x(j-1))*(rho(j)+rho(j-1)));
         end
      elseif mmpde==5
         for j=2:(jmax-1)
            Resx(j) = tau*xt(j)-0.5*dxi*dxi*((x(j+1)-x(j))*(rho(j+1)+rho(j)) ...
                 -(x(j)-x(j-1))*(rho(j)+rho(j-1)));
         end
      else % for modified MMPDE5 (default)
         for j=2:(jmax-1)
            Resx(j) = tau*xt(j)-0.5*dxi*dxi/rho(j)*((x(j+1)-x(j))*(rho(j+1)+rho(j)) ...
                 -(x(j)-x(j-1))*(rho(j)+rho(j-1)));
         end
      end
              
      % for mesh pde: boundary points (fixed)
   
      Resx(1) = xt(1)-0.0;
      Resx(jmax) = xt(jmax)-0.0;
   end
   
   res = [reshape(Res,nn,1)
          Resx];

% end of movefem_res

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

function [SPDY,SPDYP] = jpattern(npde,jmax,mmpde,alpha_type)
%
% this function defines the sparse pattern of the jacobian matrix.
%

   SPDY = sparse((npde+1)*jmax, (npde+1)*jmax);
   SPDYP = sparse((npde+1)*jmax, (npde+1)*jmax);
   nn = npde*jmax;

   for j=1:jmax
      i=(npde*(j-1)+1):(npde*j);
      j1 = max(j-1,1);
      j2 = min(j+1,jmax);	  
      SPDY(i,(npde*(j1-1)+1):(npde*j2)) = ones(npde,(j2-j1+1)*npde);
      SPDY(i,(nn+j1):(nn+j2)) = ones(npde,j2-j1+1);
   end
   SPDYP(1:nn,1:end) = SPDY(1:nn,1:end);

   if alpha_type==3  % for alpha = constant
      for j=1:jmax
	     j1 = max(j-6,1);
	     j2 = min(j+6,jmax);
	     SPDY(nn+j,(npde*(j1-1)+1):(npde*j2)) = ones(1,(j2-j1+1)*npde);
	     SPDY(nn+j,(nn+j1):(nn+j2)) = ones(1,j2-j1+1);
      end
   else % for alpha based on solution derivatives   
      SPDY((nn+1):end,:) = ones(jmax,nn+jmax);
   end

   e = ones(jmax,1);
   if (mmpde==4) | (mmpde==6) 
      SPDYP((nn+1):end,(nn+1):end) = spdiags([e e e],-1:1,jmax,jmax);
   else
      SPDYP((nn+1):end,(nn+1):end) = spdiags(e,0,jmax,jmax);
   end
   
% end of jpattern

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

function tau = TauFun_default(npde,jmax,t,x,v)

   tau = 0.01;


   
