%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[]=Program
%% the code is a copy of xie.m
%% example 1: blood vessel extension without proliferation

        global XZA   XZB   XZD  XZbeta   vis  Elast
        global XZUSE_ABS
        global ZiterMax

clf

% set parameters:
        ZiterMax = 1;  %% most iteration allowed in each time step. If =1, then no iteration
        XZUSE_ABS  = true;  %% choose to use |1+u_x| or not, its value=true OR false
        XZbeta  = 0.01;      %% beta
        vis  = 1.d-4;         %% \mu, viscosity
        Elast  = 1.d0;       %% Young's modulus for elasticity, always=1

	nx = 200;		%% spatial discretization points on [0,1]
        dt= 1.d-3;     	 	 %% time step 
        time_end = 7;  	%% final time
        dt_output =1d0; 		 %% time interval to output results


% associated parameters:
        hx = 20.d0/nx;
        x=linspace(0,20,nx+1); 
        out_freq = dt_output/dt;

% initial values:
        u=zeros(nx+1,1);
        uxp1=ones(nx+1,1);
        unew=u;
        uxp1new=uxp1;

	strg(1,:) = ' u   ';
	strg(2,:) = '1+u_x';
	%strg(3,:) = 'exact';


% start:
        time=0;      
        iter = 0;
        while time <= time_end
             [unew, uxp1new]=ElasticGrowth(u, uxp1, x, time, dt, hx, nx);
             u=unew; uxp1=uxp1new;
             iter=iter+1;
             time=time+dt;
             if mod(iter, out_freq) ==0 
                err_estimate(u, uxp1, time, dt, x, hx, nx);
                %clf; 
                hold on;
                %Han(1)=plot(x,u, 'r-','linewidth', 2);
                plot(x,u, 'r-','linewidth', 2);
                set(gca,'FontSize',18); %% set the font of tick label
                xlabel('x, Lagrangian coordinate', 'fontsize', 20);
                ylabel('Displacement (X100 \mu m)', 'fontsize', 20);
                %Han(2)=plot(x,uxp1, 'g--', 'linewidth', 2);
                %Han(3)=plot(x,uexact(time,x), 'b--', 'linewidth', 2);
	 	%IP=legend(Han,strg,3);
		%set(IP,'FontSize', 20)
                %title(sprintf('Time = %10.6E', time),'Fontsize', 20 );
                disp(sprintf('%12.5e  ', time, u(nx+1)));
                %axis([0,1,-1,1]);
                axis([0,20,-20,2]);
                pause(0.1);
             end
        end
        axis([0,20,-20,2]);
        print('-deps', 'regress4.eps')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      function[]= err_estimate(u, uxp1, time, dt, x, hx, nx);
 
      errmax = 0.d0;
      errL2 = 0.d0;
      errH1 = 0.d0;

      errmax = max(abs(u'-uexact(time,x)));
      errL2 = sqrt( sum( (u'-uexact(time,x)).^2)*hx);
      errH1 = sqrt( sum( (u'-uexact(time,x)).^2)*hx ...
                  + sum( (uxp1'-1-udexact(time, x)).^2)*hx ); 

      disp('   time          dt              Max Err      L2 Err        H1 Err');
      %disp(sprintf('%12.5e, %12.5e, %12.5e, %12.5e, %12.5e', time, dt, errmax, errL2, errH1));
      disp(sprintf('%12.5e  ', time, dt, errmax, errL2, errH1));
      disp(sprintf('\n  '));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      function[us]=crout(AD,AU,AL,b,N);

      u=zeros(N,1);
      z=zeros(N,1);
%step 4:
      t=AD(1);
      u(1) = AU(1)/AD(1);
      z(1) = b(1)/t;

%step 5:
      for i=2:N-1
          t=AD(i)-AL(i)*u(i-1);
          u(i) = AU(i)/t;
          z(i) = (b(i)-AL(i)*z(i-1))/t;
      end

%step 6:
      t=AD(N)-AL(N)*u(N-1);
      z(N)=( b(N) - AL(N)*z(N-1) )/t ;

      x(N)=z(N) ;
      for i=N-1:-1:1
          x(i) = z(i) - u(i)*x(i+1);
      end

      us=x(1:N);
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      function[unew, uxp1new]=ElasticGrowth(u, uxp1, x, time, dt, hx, nx);
        global XZA   XZB   XZD  XZbeta   vis  Elast
        global XZUSE_ABS
        global ZiterMax


      AL=zeros(2,2);
      RAL=zeros(2,2);
      vgp1D(1,1) = 1.d0;
      vgp1D(1,2) = 0.5d0;
      vgp1D(1,3) = 0.d0;
      vgp1D(2,1) = 0.d0;
      vgp1D(2,2) = 0.5d0;
      vgp1D(2,3) = 1.d0;
      gw1D(1) = 1.d0/6;
      gw1D(2) = 2.d0/3;
      gw1D(3) = 1.d0/6;

      gp(1) = -1/hx;
      gp(2) = 1/hx;

      lambda = uxp1;
      uold   = u;
%%!! iteration start:
      err = 1.d0;
      iter = 0;
      while err > 1.d-10 & iter < ZiterMax
         if XZUSE_ABS
            lamiter = abs(lambda);
         else
            lamiter = lambda;
         end
         rho = growth(time, x);
         for k=1:nx
            XBD(k) = 0.d0;
            XBU(k) = 0.d0;
            XBL(k) = 0.d0;
            XF(k) = 0.d0;
         end
         for k=1:nx
            rgp1D(1) = hx*(k-1);
            rgp1D(2) = hx*( 2*k-1 )/2;
            rgp1D(3) = hx*(k);
          
            for i=1:3
               lamiterl(i) = vgp1D(1,i)*lamiter(k) + vgp1D(2,i)*lamiter(k+1);
               rhol(i) = vgp1D(1,i)*rho(k) + vgp1D(2,i)*rho(k+1);
            end

            %%!! local matrix
            %%!! local matrix
            %%!! local matrix
            for i=1:2
                for j=1:2
                  AL(i,j) = 0.d0;
                  RAL(i,j) = 0.d0;
                  for l=1:3
                     tmp1 = XZbeta/dt*vgp1D(i,l)*vgp1D(j,l)*lamiterl(l);

                     tmp2 = (vis/dt+ Elast)*gp(i)*gp(j);
                     AL(i,j) = AL(i,j) + gw1D(l)*(tmp1 + tmp2);

                     tmp2 = (vis/dt)*gp(i)*gp(j);
                     RAL(i,j) = RAL(i,j) + gw1D(l)*(tmp1 + tmp2);
                  end
                  AL(i,j) = AL(i,j)*hx;
                  RAL(i,j) = RAL(i,j)*hx;
               end
            end

            %%!! local load
            for i=1:2
               b(i) = 0.d0;
               for j=1:2
                  b(i) = b(i) + RAL(i,j)*uold(k-1+j);
               end
               for l=1:3
                  b(i) = b(i) + gw1D(l)*rhol(l)*gp(i)*hx;
               end
            end

            %%!! assemble matrix and load:
            if k == 1
               XBD(1) = AL(2,2);
               XF(k) = XF(k) + b(2);
            else
               XBD(k-1) = XBD(k-1) + AL(1,1);
               XBD(k) = XBD(k) + AL(2,2);
               XBL(k) = AL(2,1);
               XBU(k-1) = AL(1,2);
               XF(k-1) = XF(k-1) + b(1);
               XF(k) = XF(k) + b(2);
            end
         end

         %%!! need to retrieve 2D protrusion force:
         XF(nx) = XF(nx) + ProtruFunc(time);

         %%!! solve:
         us = zeros(nx+1,1);
         KCF = zeros(1,nx);
         KCF=crout(XBD,XBU,XBL,XF,nx);
         us(2:nx+1)=KCF(1:nx);
         us(1) = 0.d0;
         for i=2:nx
            lambda(i) = ( us(i+1)- us(i-1) ) /2/hx+ 1.d0;
         end
         lambda(1) = ( -3*us(1) + 4*us(2)- us(3) ) /2/hx+ 1.d0;
         lambda(nx+1) =( us(nx-1) -4*us(nx)+3*us(nx+1) ) /2/hx+ 1.d0;

         %%!! iterative error:
         err = max(abs(lambda-lamiter));
         %disp(sprintf('time=%5.3f, err= %5.4e\n',time, err));

         %%!! update sprout information:
         uxp1 = lambda;
         u = us;
         iter = iter + 1;
      end

      uxp1new = lambda;
      unew = us;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Do not modify this 'udexact'  %%%%%%%%%%%%%%%%%%%%%%%
      function[ud]=udexact(t, x);
      global XZA   XZB   XZD  XZbeta   vis  Elast

      ud=   growth(t,x)-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      function[u]=uexact(t, x);
      global XZA   XZB   XZD  XZbeta   vis  Elast

      u= XZA/3*x.^3 + (XZB-1)*x;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% g function
      function[g]= ProtruFunc(t);
      global XZA   XZB   XZD  XZbeta   vis  Elast

      g = 0.d0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% f function
      function[f]= growth(t, x);
      global XZA   XZB   XZD  XZbeta   vis  Elast

        XZA  = 1;         %% some parameter, such as XZA x^2
        XZB  = 1;         %% some parameter, such as XZA x^2
        XZD  = 0;         %% some parameter, such as XZA x^2
      %[n1,n2]=size(x);
      %f=zeros(n1,n2);
      %for i=1:n2
      %    f(1,i)= max(0,1-t/7);
      %end
      f= max(0,1-t/7-t*x/140);
      f= (1-t/7)*cos(x/20);

      f = f -1.d0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
