function fem = fpc(Reynolds,Blockage,Tmax,pert,PLOTS,MOVIE)

% function fpc(Reynolds,Blockage,Tmax,pert,PLOTS,MOVIE)
%
% Solves flow past a cylinder for problem definition given in 
% fem_template_fpc. 
%
% First solves the steady state problem, then uses the
% solution to initialize the time dependent simulation 
% on the time interval linspace(0,Tmax,10).
%
% take as inputs :
% Reynolds number
% Blockage ratio
% Tmax = max time in timedependent solve
% pert = boundary perturbation to create asymmetry.
% PLOTS and MOVIES are flags, either (0/off or 1/on) 
% on the time interval [0:Tmax], recording 10 intermediate solutions.
% if some inputs are missing, set to defaults (10,0.5,0.5,0.05,0,0)
%
% The entire fem structure as saved as fpc_R(Re#)_B(B).mat
% which can then be loaded with matlab and viewed.
%
% The equation that fem_template_fpc solves is
% NS) rho u' - div(eta (grad(u)+grad(u)^t)) + rho (u.grad u) + grad(p) = F
% with
% rho = Re
% eta = 1
% F = 0 (could include bouyancy here)
%
% fem_template_fpc uses as boundary conditions:
% parabolic inflow, Umax*4*s*(1-s) , (max flow = Umax)
% with Umax = 1/(1-B^2/3), forcing Ubar = 1
% p=0 at outflow
% and no slip otherwise 
%

close all;

load template_fpc 

% Do some input checking
% if some inputs are missing, set to defaults (10,.5,.5,0.05)
switch nargin
   case 0
      Reynolds=10;Blockage=0.5;Tmax=1/Blockage;pert=0.05;PLOTS=0;MOVIE=0;
   case 1
      Blockage=0.5;Tmax=1/Blockage;pert=0.05;PLOTS=0;MOVIE=0;
   case 2
      Tmax=1/Blockage;pert=0.05;PLOTS=0;MOVIE=0;
   case 3
      pert=0.05;PLOTS=0;MOVIE=0;
   case 4
      PLOTS=0;MOVIE=0;
  otherwise 
      MOVIE=0;
end

% check if input is physical
if (Reynolds <0 | Blockage>=1 | Blockage<=0 | Tmax <= 0)
    error('You entered unphysical parameters')
end

% equation parameters
Re = Reynolds;
% Geometry  paramters
D = 1;                  % Diameter of pipe
d = Blockage*D;         % diameter of cylinder
B = Blockage;           % blockage ratio (d/D)
Umax = 1/(1-B^2/3);     % reciprocal of Ubar with max inflow = 1

% Constant used only in solve step
flowmax = 1/(1-B);
fem.flowmax = flowmax;

% Overwrite the constants.  No others are used by fem_template_fpc 
fem.const = {'Re' Re 'D' D 'd' d 'B' B 'Umax' Umax};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start to reset selected parts of template
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% set geometry to -3<x<6  0<y<D+pert - circle w/ c (0,D/2) and r = d/2 
% Use a smaller asymmetry as blockage and flowrate increase?
fem.geom = rect2(-3,6,0,D+pert) - circ2(0,D/2,d/2);

% Initialize the mesh
fem.mesh = meshinit(fem,'hmax',0.1,'hgrad',1.1,'hcurve',0.1);

% Extend the mesh
fem.xmesh = meshextend(fem);

if PLOTS == 1
   % Info about and plot the current mesh
   figure;
   meshinfo(fem.mesh)
   meshplot(fem); axis equal;
end

% Compute steady state soln
% Solve the equations and plot the result.
fem.sol = femnlin(fem,'report','on','stop','on','Ntol',1e-3*flowmax);
%fem.sol = femnlin(fem,solncomp',{'u'},'report','on','stop','on','Ntol',1e-4);

% save for restart
fem0 = fem;

if PLOTS == 1
   figure;
   % Plot initial solution
   postplot(fem,'tridata','sqrt(u^2+v^2)',...
             'tribar','on',...
             'contdata','p',...
             'contbar','on',...
             'flowdata',{ 'u' 'v'},...
             'title', 't = 0 Surface: Velocity (U) Contour: pressure (p)  Flow: [x velocity (u),y velocity (v)]  ',...
             'solnum', 1,...
             'axisequal','on');

   figure;
   % Plot final time solution
   postplot(fem,'tridata','sqrt(u^2+v^2)',...
             'tribar','on',...
             'contdata','p',...
             'contbar','on',...
             'flowdata',{ 'u' 'v'},...
             'title', 'Surface: Velocity (U) Contour: pressure (p)  Flow: [x velocity (u),y velocity (v)]  ',...
             'axisequal','on');
end

% move cylinder to centerline 

fem.geom = rect2(-3,6,0,D) - circ2(0,D/2,d/2);

% (Re)Initialize the mesh for the new geometry
fem.mesh = meshinit(fem,'hmax',0.1,'hgrad',1.2,'hcurve',0.1);

% Extend the mesh
fem.xmesh=meshextend(fem);

if PLOTS == 1
   figure;
   % plot mesh and give info
   meshinfo(fem.mesh)
   meshplot(fem); axis equal;
end

% Map initial solution to current xmesh
% clears old fem0.sol.u
fem0.sol.u(:,1:end-1)=[];

umap={'p',1,1,1,1;'u',1,1,1,1;'v',1,1,1,1};
init=asseminit(fem,'init', {fem0,umap});

% Restart and compute a time-dependent simulation
% to look for Karman vortex trails.

% Solve dynamic problem
fem.sol = femtime(fem,'tlist',linspace(0,Tmax,10),...
       'ode','fldaspk',...
       'atol',{'u' 1e-3*flowmax 'v' 1e-3*flowmax 'p' inf },...
       'stop' , 'on',...
       'init',init,...
       'timeind','D',...
       'report','on');
% stop on returns partial solution of timestepping fails
% init on uses initial condition (built by assemint) 
% timeind D, mass matrix is constant (for speed)
% see ref manual pg 5-134 for full description of femtime command

% Save current fem structure for restart purposes
fem0=fem;

% Establish a distinct filename
fname = sprintf('fpc_R%d_B%d_t%3.2f',Re,B*100,fem.sol.tlist(end));

% Generate some graphics if the flags are set

if PLOTS == 1

   figure;
   % Plot initial solution
   postplot(fem,'tridata','sqrt(u^2+v^2)',...
             'tribar','on',...
             'contdata','p',...
             'contbar','on',...
             'flowdata',{ 'u' 'v'},...
				 'title', 't = 0 Surface: Velocity (U) Contour: pressure (p) Flow: [x velocity (u),y velocity (v)]  ',...
             'solnum', 1,...
             'axisequal','on');
   figure;
   % Plot final time solution
   postplot(fem,'tridata','sqrt(u^2+v^2)',...
             'tribar','on',...
             'contdata','p',...
             'contbar','on',...
             'flowdata',{ 'u' 'v'},...
             'title', 'Surface: Velocity (U) Contour: pressure (p)  Flow: [x velocity (u),y velocity (v)]  ',...
             'axisequal','on');
end

if MOVIE == 1;
   figure;
   fem.movie = postmovie(fem,'tridata','sqrt(u^2+v^2)',...
                     'tribar','on',...
                     'contdata','p',...
                     'contbar','on',...
                     'flowdata',{'u' 'v'},...
                     'title', 'Surface: Velocity (U) Contour: pressure (p)  Flow: [x velocity (u),y velocity (v)]  ',...
                     'axisequal','on',...
                     'filename',strcat(fname,'.mpg'),...
                     'filetype','mpeg');
end

save(strcat(fname,'.mat'),'fem'); 

