function fphc(Reynolds,Prandtl,Gamma,Blockage,Tmax,Tstep,pert,PLOTS,MOVIE)

% function fphc(Reynolds,Prandtl,Gamma,Blockage,Tmax,Tstep,pert,PLOTS,MOVIE)
% fem_template_fphc is the multiphysics setup
% for flow past a heated cylinder with
% Reynolds number
% Prandtl number 
% Gamma is the coupling coeff of viscosity to temperature
% Blockage ratio
% Tmax = max time in time dependent solve
% Tstep = # of time step in solve
% pert = boundary perturbation to create asymmetry.
% PLOTS and MOVIE are flags (0 | 1) to activate 
%
% The equation that fem_template_fphc uses is
% NS) rho u' - div(eta (grad(u)+grad(u)^t)) + rho (u.grad u) + grad(p) = F
% CD) Mass u' - div( c grad(T) ) + (beta. grad(T)) = f
%
% with
%
% rho = Re
% eta = 1-Gamma*T
% F = 0 (could include bouyancy here)
% Mass = 1               (denoted by d_a in femlab)
% c = 1/(Re*Pr)
% beta = [u v]
% f = 0
%
% 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
% no slip otherwise 
% Temp = T0 on outer boundary
% Temp = T1 in cylinder

close all;

load template_fphc 

% Do some input checking
% if some are missing, set them to the defaults
% Reynolds=10;Prandtl=1;Gamma=0;Blockage=0.5;Tmax=1;Tstep=3;pert=0.05;
switch nargin
   case 0
     Reynolds=10;Prandtl=1;Gamma=0;Blockage=0.5;Tmax=1;Tstep=3;pert=0.05;PLOTS=0;MOVIE=0; 
   case 1
     Prandtl=1;Gamma=0;Blockage=0.5;Tmax=1;Tstep=3;pert=0.05;PLOTS=0;MOVIE=0;
   case 2
     Gamma=0;Blockage=0.5;Tmax=1;Tstep=3;pert=0.05;PLOTS=0;MOVIE=0;
   case 3
     Blockage=0.5;Tmax=1;Tstep=3;pert=0.05;;PLOTS=0;MOVIE=0;
   case 4
     Tmax=1;Tstep=3;pert=0.05;PLOTS=0;MOVIE=0;
   case 5
     Tstep=3;pert=0.05;PLOTS=0;MOVIE=0;
   case 6 
     pert=0.05;PLOTS=0;MOVIE=0;
   case 7 
     PLOTS=0;MOVIE=0;
   otherwise 
     MOVIE=0;
end

if (Reynolds <0 | Prandtl < 0 | Blockage>=1 | Blockage<=0 | Tmax <= 0)
    error('You entered unphysical parameters')
end

% equation parameters
Re = Reynolds;
Pr = Prandtl;
gamma0 = Gamma;
T0 = 0;
T1 = 1;

% 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

% not used in fem.const
flowmax = 1/(1-B);
fem.flowmax = flowmax;

% Overwrite the constants.  No others are used by fem_template_fphc 
fem.const = {'Re' Re 'Pr' Pr 'gamma0' gamma0 'T0' T0 'T1' T1 'D' D 'd' d 'B' B 'Umax' Umax};

% 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 increases?
fem.geom = rect2(-3,6,0,D+pert) - circ2(0,D/2,d/2);

% Initialize the mesh
fem.mesh = meshinit(fem,'hmax',0.2,'hgrad',1.2,'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);

% save for restart
fem0 = fem;

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

   figure;
   % Plot final time solution
   postplot(fem,'tridata','T',...
             'tribar','on',...
             'contdata','p',...
             'contbar','on',...
             'flowdata',{ 'u' 'v'},...
             'title', 'Surface: Temp (T) 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);

% Multiphysics
% I don't know if this is needed but it doesn't seem to hurt
fem=multiphysics(fem);

% (Re)Initialize the mesh for the new geometry
fem.mesh = meshinit(fem,'hmax',0.2,'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
fem0.sol.u(:,1:end-1)=[]; %clears old fem0.sol.u
umap={'T',1,1,1,1;'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 see the Karman vortex trails.

% Solve dynamic problem
% 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

tlist = linspace(0,Tmax,Tstep);
fem.sol = femtime(fem,'tlist',tlist,...
       'ode','fldaspk',...
       'atol',{'u' 1e-3*flowmax 'v' 1e-3*flowmax 'p' inf 'T' 1e-3},...
       'stop' , 'on',...
       'init',init,...
       'timeind','D',...
       'report','on');

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

if PLOTS == 1

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

if MOVIE == 1;
   figure;
   Mov = postmovie(fem,'tridata','T',...
                     'tribar','on',...
                     'contdata','p',...
                     'contbar','on',...
                     'flowdata',{'u' 'v'},...
                     'title', 'Surface: Temp (T) Contour: pressure (p)  Flow: [x velocity (u),y velocity (v)]  ',...
                     'axisequal','on');
end

fname = sprintf('fphc_R%d_G%d_B%d_t%f.mat',Re,Gamma,B*100,Tmax);
save(fname,'fem');

