% 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;

% Geometry
g1=rect2('2','1','base','corner','pos',{'-1','-1'},'rot','0');

% Analyzed geometry
clear s
s.objs={g1};
s.name={'R1'};
s.tags={'g1'};

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

% Initialize mesh
fem.mesh=meshinit(fem, ...
    'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'FlNavierStokes';
appl.gporder = {4,2};
appl.cporder = {2,1};
appl.assignsuffix = '_ns';
clear prop
prop.analysis='static';
appl.prop = prop;
clear bnd
bnd.type = {'open','stress','walltype'};
bnd.stresstype = {'ntotnflow','ntotnflow','ntotnflow'};
bnd.walltype = {'noslip','noslip','lwall'};
bnd.opentype = {'ntotstress','ntotstress','novisc'};
bnd.uwall = {0,0,'Vs'};
bnd.ind = [2,2,3,2];
appl.bnd = bnd;
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;

% Subdomain settings
clear equ
equ.ind = [1];
equ.dim = {'u','v','p'};


% Subdomain settings
clear equ
equ.ind = [1];
equ.dim = {'u','v','p'};

% Subdomain expressions
equ.expr = {'Vs','-0.5+(x>0)',...
    'r','sqrt(x^2+y^2)', ... radius from top left corner
    'th','asin(x/r)', ... angle from vertical
    'pc','-4*cos(th)/(pi*r)', ... pressure solution
    'prc','-4*cos(th)/pi', ... pressure x radius solution
    'vxc','2*(th-sin(th)*cos(th))/pi', ... x-velocity
    'vyc','2*cos(th)^2/pi'}; %y-velocity

fem.equ = equ;
% 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
fem=adaption(fem, ...
    'blocksize',1000, ...
    'solcomp',{'u','p','v'}, ...
    'outcomp',{'u','p','v'}, ...
    'hnlin','on', ...
    'solver','stationary', ...
    'l2scale',[1], ...
    'l2staborder',[2], ...
    'eigselect',[1], ...
    'maxt',10000000, ...
    'ngen',2, ...
    'resorder',[0], ...
    'rmethod','meshinit', ...
    'tppar',1.7, ...
    'uscale','none', ...
    'nonlin','off',...
    'geomnum',1);

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

%% Plot solution
zm=0.075;
figure(1); clf;
hold on;
Hm=meshplot(fem,...
    'edgecolor',[1,1,1]*3/4)
Hf=postplot(fem, ...
    'flowdata',{'u','v'}, ...
    'flowdens','uniform', ...
    'flowcolor',[0.0,1.0,0.0], ...
    'flowlines',20,...
    'geom','off');
Ha=postplot(fem, ...
    'arrowdata',{'u','v'}, ...
    'arrowxspacing',21, ...
    'arrowyspacing',11, ...
    'arrowtype','arrow', ...
    'arrowstyle','proportional', ...
    'arrowcolor',[1,0.0,1]*3/4, ...
    'geom','off');
set(Hf,'linewidth',4);
hold on;plot([1,1,-1,-1,1]*zm,[0,-1,-1,0,0]*zm,'k','linewidth',2)

axis equal; 
axis tight; 
A1=gca;

A2=axes('position',[0.0,0.1,0.6,0.3]);
postplot(fem, ...
    'tridata',{'u','cont','internal'}, ...
    'trimap','jet(1024)', ......
    'tribar','off',...
    'geom','off'); 
hold on;
meshplot(fem,...
    'edgecolor',[1,1,1]*3/4);
hold on;
plot([1,1,-1,-1,1]*zm,[0,-1,-1,0,0]*zm,'k','linewidth',2)
set(gca,'xtick',[],'ytick',[])
axis equal; 
axis([-1,1,-1,0]*zm); 

