% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)
% Some geometry objects are stored in a separate file.
% The name of this file is given by the variable 'flbinaryfile'.

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;

% initalize time counter...
tic
% Geometry 1x1 rectangle
g1=rect2('1','1','base','corner','pos',{'0','-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: free mesh with quadrangles. Versions before 3.4 are not
% ble to do that. Just comment out the methodsub line to have a free mesh
% with triangles
fem.mesh=meshinit(fem, ...
    'hauto',4,...    
    'hmaxvtx',[2,0.000005], ...
    'methodsub','quad');

% Application mode 1: Navier Stokes
clear appl
appl.mode.class = 'FlNavierStokes';
% Lagrange P2 P1 element (important for pressure) 
appl.shape = {'shlag(2,''u'')','shlag(2,''v'')','shlag(1,''p'')'};
appl.gporder = {4,2};
appl.cporder = {2,1};
appl.assignsuffix = '_ns';
clear prop
prop.analysis='static';
appl.prop = prop;
clear pnt
pnt.pnton = {0,1};
pnt.p0 = {0,'pc'}; % point pressure BC
pnt.ind = [2,1,2,2]; %fix P  at 3 corners
appl.pnt = pnt;
clear bnd
% 1: symmetry; 2: analytical flow field; 3: uniform velocity
bnd.type = {'sym','walltype','walltype'};
bnd.walltype = {'noslip','lwall','lwall'};
bnd.uwall = {0,'vxc',1};
bnd.vwall = {0,'vyc',0};
bnd.ind = [1,2,3,2];
appl.bnd = bnd;
clear equ1
equ.cporder = {{1;1;2}};
equ.gporder = {{1;1;2}};
equ.rho = 0;
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 expressions
equ.expr = {'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;
t1=toc; %stop time counter to register setup time
tic; % restart time counter
% Multiphysics
fem=multiphysics(fem);
% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femlin(fem, ...
    'solcomp',{'u','p','v'}, ...
    'outcomp',{'u','p','v'});

% % Solve problem (AMR)
% fem=adaption(fem, ...
%     'solcomp',{'u','p','v'}, ...
%     'outcomp',{'u','p','v'}, ...
%     'solver','stationary', ...
%     'l2scale',[1], ...
%     'l2staborder',[2], ...
%     'eigselect',[1], ...
%     'maxt',10000000, ...
%     'ngen',1, ...
%     'resorder',[0], ...
%     'rmethod','meshinit', ...
%     'tpfun','fltpworst', ...
%     'tppar',0.5, ...
%     'geomnum',1);

t2=toc; % Record time to solution

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

%% Calculate global error
% Integrate
I1=postint(fem,'((p/pc-1)^2)', ...
    'dl',[1]);
disp(sprintf('Time to solution: %g+%g=%gs; %g d.o.f.',t1,t2,t1+t2,flngdof(fem)));
disp(sprintf('Global error: %g on pressure',I1));
     
%% Plot solution
figure(1); clf;
subplot(221);
meshplot(fem);
axis equal
axis([0,1,-1,0])

subplot(222);
Hg=plot([0,1,1,0,0],-[0,0,1,1,0],'k','linewidth',1);
hold on;
Hp=postplot(fem, ...
    'contdata',{'p','cont','internal','unit','Pa'}, ...
    'contlevels',[[-20:2:0]], ...
    'contlabel','off', ...
    'contmap','jet(1024)','geom','off');
set(Hp,'linewidth',2);
Hpc=postplot(fem, ...
    'contdata',{'pc','cont','internal','unit','Pa'}, ...
    'contlevels',[[-20:2:0]], ...
    'contlabel','off', ...
    'contbar','off',...
    'contmap','jet(1024)','geom','off');
set(Hpc,'linewidth',0.5);
Hac=postplot(fem, ...
    'arrowdata',{'u','v'}, ...
    'arrowxspacing',7, ...
    'arrowyspacing',7, ...
    'arrowtype','arrow', ...
    'arrowstyle','proportional', ...
    'arrowcolor',[1.0,0.5,1.0], ...
    'geom','off');
set(Hac,'linewidth',2)
Hac=postplot(fem, ...
    'arrowdata',{'vxc','vyc'}, ...
    'arrowxspacing',7, ...
    'arrowyspacing',7, ...
    'arrowtype','arrow', ...
    'arrowstyle','proportional', ...
    'arrowcolor',[0.5,0.0,0.5], ...
    'geom','off');
set(Hac,'linewidth',0.5)
axis equal
axis([0,1,-1,0])

%% Plot in cross-section or along domain
subplot(212);
postcrossplot(fem,1,[1], ...
    'lindata','-p*r*pi/4', ...
    'cont','internal', ...
    'linxdata',{'-y','unit','m'}, ...
    'refine','auto');
xlabel('Depth','fontsize',18);
ylabel('-p r \pi /4','fontsize',18);
axis([0,1,0.5,1.5])
grid off;
box on;

title(sprintf('Time to solution: %gs; int((p/p_c-1)^2)=%g; %g dof',t2,I1,flngdof(fem)));
