% Re du/dt +dp/dx-(1-alpha) d2u/dx2 - dtau/dx=F, (2)
% dp/dt+beta du/dx = beta R,        		     (1)
% dtau/dt+tau1De-2(tau+alpha/De)du/dx = 0        (3)

% u is fluid velocity field,
% p is fluid pressure, and
% tau is the stress tensor.
% F is the body force.

% The boundary conditions are no-flux at all boundaries,
% no-slip for the velocity
%
clear fem1 fem2

L=6.0; % Domain length
beta=0.3; % inverse compressibility
Re=0.1; % Reynolds number
De=0.1;; % Deborah number
alpha=0.9; % Non-Newtonian fraction of viscosity
dt=0.05; % time step
t=0; % global time variable
xleft=-L/2; % current position of the left end point
xright=L/2; % current position of the right end point

k=0; % Initialize counter
kend=5; % number of steps


% working structure
fem1.sdim={'x'};
fem1.form = 'coefficient';  
fem1.dim = {'u','p','tau'}; 
fem1.shape=2;

% The main loop starts

while k<kend
    
k=k+1;


fem1.const = {'LL',L,'Beta',beta,'RRe',Re,'DDe',De,'al',alpha};
fem1.expr = {'F' '-0.2*x' 'R' '0'};
% domain geometry
r=solid1([xleft xright]);
% Define objects
objs={r};
names={'r'};
s.objs=objs;
s.name=names;
% 
% Curve objects
%
objs={};
names={};
c.objs=objs;
c.name=names;

objs={};
names={};
p.objs=objs;
p.name=names;

drawstruct=struct('s',s,'c',c,'p',p);

fem1.draw=drawstruct;
fem1.geom=geomcsg(fem1);

fem1.mesh=meshinit(fem1);

fem1.equ.a = {{0 0 0; ...
            0 0 0; ...
            0 0 1 ...
        }};
fem1.equ.c = {{'2*(1-al)' 0 0}};
fem1.equ.be = {{ ...
0 0 0; ...
'Beta' 0 0; ...
'-2*(tau*DDe+al)' 0 0 ...
}};  

fem1.equ.al = {{ ...
 0 -1 1; ...
 0 0 0; ...
 0 0 0 ...
}};    
fem1.equ.f={{...
             'F' ...
             'Beta*R' ...
             '0' ...
         }};
 
%fem1.bnd.q = {...
%{'dnx' 'dny' '0'; ...
% 'dnx' 'dny' '0'; ...
%  '0'   '0'  '0'} ...
%};
%fem1.bnd.g = {{'0' ...
%              '0' ...
%               '0'}};
   
fem1.equ.da = {{'RRe' 1 'DDe'}};

if k==1
fem1.equ.init={{'0.0' '0.0' '0'}};
else 
fem1.equ.init={{'femval_stat1(x,''u'',''fem2'')' ...
        'femval_stat1(x,''p'',''fem2'')' ...
    'femval_stat1(x,''tau'',''fem2'')' ...
}};
end

fem1.xmesh = meshextend(fem1);

% fem.sol = femlin(fem);
fem1.sol = femtime(fem1,'tlist',[0:dt/5:dt],'report','off');

%postplot(fem,'lindata','u','linbar','on');
%postplot(fem,'liny','u','linstyle','r'); hold on;
%postplot(fem,'liny','p','linstyle','b'); hold on;
%postplot(fem,'liny','tau','linstyle','g'); hold off;

% compute the velocities of the end points
% using postinterp
uend=postinterp(fem1,'u',[xleft xright]);
% find the end points shift xleft xright
xleft = xleft+dt*uend(1);
xright=xright+dt*uend(2);
uaver=(abs(uend(1))+abs(uend(2)))/2;
% update time 
t=t+dt;

[t xright uaver]

fem2.sdim={'x'};
fem2.form = 'coefficient';  
fem2.dim = {'u','p','tau'}; 
fem2.shape=2;

fem2.const = {'LL',L,'Beta',beta,'RRe',Re,'DDe',De,'al',alpha};
fem2.expr = {'F' '-0.2*x' 'R' '0'};
% domain geometry
r=solid1([xleft xright]);
% Define objects
objs={r};
names={'r'};
s.objs=objs;
s.name=names;
% 
% Curve objects
%
objs={};
names={};
c.objs=objs;
c.name=names;

objs={};
names={};
p.objs=objs;
p.name=names;

drawstruct=struct('s',s,'c',c,'p',p);

fem2.draw=drawstruct;
fem2.geom=geomcsg(fem2);

fem2.mesh=meshinit(fem2);

fem2.equ.a = {{0 0 0; ...
            0 0 0; ...
            0 0 1 ...
        }};
fem2.equ.c = {{'2*(1-al)' 0 0}};
fem2.equ.be = {{ ...
0 0 0; ...
'Beta' 0 0; ...
'-2*(tau*DDe+al)' 0 0 ...
}};  

fem2.equ.al = {{ ...
 0 -1 1; ...
 0 0 0; ...
 0 0 0 ...
}};    
fem2.equ.f={{...
             'F' ...
             'Beta*R' ...
             '0' ...
         }};
 
%fem2.bnd.q = {...
%{'dnx' 'dny' '0'; ...
% 'dnx' 'dny' '0'; ...
%  '0'   '0'  '0'} ...
%};
%fem2.bnd.g = {{'0' ...
%              '0' ...
%               '0'}};
   
fem2.equ.da = {{'RRe' 1 'DDe'}};

fem2.equ.init={{'femval_stat1(x,''u'',''fem1'')' ...
        'femval_stat1(x,''p'',''fem1'')' ...
    'femval_stat1(x,''tau'',''fem1'')' ...
}};


fem2.xmesh = meshextend(fem2);

% fem.sol = femlin(fem);
fem2.sol = femtime(fem2,'tlist',[0:dt/5:dt],'report','off');

%postplot(fem,'lindata','u','linbar','on');
%postplot(fem,'liny','u','linstyle','r'); hold on;
%postplot(fem,'liny','p','linstyle','b'); hold on;
%postplot(fem,'liny','tau','linstyle','g'); hold off;

% compute the velocities of the end points
% using postinterp
uend=postinterp(fem2,'u',[xleft xright]);
% find the end points shift xleft xright
xleft = xleft+dt*uend(1);
xright=xright+dt*uend(2);
uaver=(abs(uend(1))+abs(uend(2)))/2;
% update time 
t=t+dt;

[t xright uaver]

end % while loop closes

postplot(fem2,'liny','u','linstyle','r'); hold on;
postplot(fem2,'liny','p','linstyle','b'); hold on;
postplot(fem2,'liny','tau','linstyle','g'); hold off;