% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.4';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 248;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2007/10/10 16:07:51 $';
fem.version = vrsn;

% Setup constants
fem.const.L=10; %length of box
fem.const.Vp=1; %velocity of top plate
fem.const.DPDx=6*fem.const.Vp; %horizontal pressure gradient

% Geometry: setup box length L, height 1, bottom left corner at [0,0]
g1=rect2(fem.const.L,'1',....
    'base','corner','pos',{'0','0'},....
    'rot','0');

% Analyzed geometry: translate into geom construct; Only object is 2D
% surface
clear s
s.objs={g1};
s.name={'R1'};
s.tags={'g1'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Create mapped quad mesh; Default is coarser than here
fem.mesh=meshmap(fem, 'hauto',2);

% (Default values are not included); e.g. viscosity=1, density=1;

% Application mode 1: Navier-Stokes
clear appl
appl.mode.class = 'FlNavierStokes';
appl.gporder = {4,2};
appl.cporder = {2,1};
appl.assignsuffix = '_ns';
% Options for application: static analysis
clear prop
prop.analysis='static';
appl.prop = prop;
% Boundary conditions
% type 1: no slip (fixed wall);
% type 2: moving/leaking wall, velocity [Vp,0];
% type 3: Impose normal stress = DPDx * x, normal flow only
clear bnd
bnd.type =       { 'walltype', 'walltype',   'stress'};
bnd.stresstype = {'totstress','totstress','ntotnflow'};
bnd.f0 =         {          0,          0,   'DPDx*x'};
bnd.walltype =   {   'noslip',    'lwall',    'noslip'};
bnd.uwall =      {          0,       'Vp',           0};
% vwall remains at default v=0 for all boundaries
% Assign each boundary type to each boundary
% boundary ID: 1,2,3,4
bnd.ind =     [3,1,2,3]; %boundary type ID
appl.bnd = bnd;
% Specifies order of element; here is default
clear equ
equ.cporder = {{1;1;2}};
equ.gporder = {{1;1;2}};
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

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

% Solve problem
% Modify default: problem is linear  use default direct solver
fem.sol=femstatic(fem, ...
    'blocksize',1000, ...
    'solcomp',{'u','p','v'}, ...
    'outcomp',{'u','p','v'}, ...
    'nonlin','off',...
    'uscale','none');

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

%% Plot solution: surface is pressure solution; arrows are velocity;
postplot(fem, ...
    'tridata',{'max(min(p-DPDx*x,1e-12),-1e-12)','cont','internal','unit',''}, ...
    'trimap','jet(1024)', ......
    'tribar','off',...
    'arrowdata',{'u','v'}, ...
    'arrowxspacing',15, ...
    'arrowyspacing',15, ...
    'arrowtype','arrow', ...
    'arrowstyle','proportional', ...
    'arrowcolor',[0.0,0.0,0.0]);
axis equal;
axis ([0,fem.const.L,0,1])
%set(gca,'position',[0,0,1000,200])
% print(1,'-dpdf','CouettePoisseuille')

%% Integrate horizontal velocity on right boundary: dl=4,edim=1 means 4th 1D
% object
I1=postint(fem,'u', 'dl',[4], 'edim',1);
% Average horizontal velocity: dl=1,edim=2 means 1st 2D object
I2=postint(fem,'u', 'dl',[1], 'edim',2)/fem.const.L;
% Show the result
disp(sprintf('Analytical flux        = %g',fem.const.Vp/2-fem.const.DPDx/12));
disp(sprintf('Flux on right boundary = %g',I1));
disp(sprintf('Average flux           = %g',I2));


