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

% Constants
fem.const.n=3; %initial stress
fem.const.L=10; % length of the box
fem.const.etamax=1e4; %lower viscosity cutoff
fem.const.etamin=1e-4; %lower viscosity cutoff

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

% Geometry: only 2D surface
% Analyzed geometry
clear s
s.objs={g1};
s.name={'R1'};
s.tags={'g1'};
fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Create mapped quad mesh: not "finer" parameters (hauto=2)
fem.mesh=meshmap(fem, 'hauto',2);

% Application mode 1: Navier Stokes
clear appl
appl.mode.class = 'FlNavierStokes';
appl.gporder = {4,2};
appl.cporder = {2,1};
appl.assignsuffix = '_ns';
% Type of analysis is "static"
clear prop
prop.analysis='static';
appl.prop = prop;
% setup boundary conditions
clear bnd
% Boundary conditions
% type 1: no slip (fixed wall);
% type 2: symmetry;
% type 3: Impose normal stress = x, normal flow only
clear bnd
bnd.type =       { 'walltype',      'sym',   'stress'};
bnd.stresstype = {'totstress','totstress','ntotnflow'};
bnd.f0 =         {          0,          0,        'x'};
bnd.ind = [3,2,1,3];
appl.bnd = bnd;
clear equ
equ.init = {{'y^2-1/4';0;'x';0;0}};
equ.eta = 'eta';
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 expressions
equ.expr = {'exy','0.5*diff(u,y)', ...
    'eta','max(min(0.5*abs(exy)^(-1+1/n),etamax),etamin)'};
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 initial case: linear (n=1) for simplicity
fem.sol=femstatic(fem, ...
    'blocksize',1000, ...
    'solcomp',{'u','p','v'}, ...
    'outcomp',{'u','p','v'}, ...
    'nonlin','on',...
    'uscale','none');
% % Solve problem
% fem.sol=femstatic(fem, ... 
%     'blocksize',1000, ...
%     'solcomp',{'u','p','v'}, ...
%     'outcomp',{'u','p','v'}, ...
%     'uscale','none');
ans
% Save current fem structure for restart purposes
fem0=fem;
femstart=fem;
% %% increase stress exponent progressive
% for i=2:ni;
%     fem.const.n=nall(i); %update stress exponent
%     disp(sprintf('Currently n=%g',fem.const.n)); %check
%     
%     % rebuild model
%     fem=multiphysics(fem); 
%     fem.xmesh=meshextend(fem);
%     % Solve: iterative solver, 
%     fem.sol=femstatic(fem, ...
%         'init',fem0.sol,... %use currently stored solution as starting guess 
%         'blocksize',1000, ...
%         'solcomp',{'u','p','v'}, ...
%         'outcomp',{'u','p','v'}, ...
%         'uscale','none');
%     fem0=fem; %store new solution 
% end
    

%% Plot solution: surface is difference between strain rate and analytical;
% arrows are velocity;
postplot(fem, ...
    'tridata',{'diff(u,y)-2*y^n','cont','internal'}, ...
    'trimap','jet(1024)', ......
    'tribar','on',...
    '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])
% print(1,'-dpdf','PowerPoiseuille')

%% Average horizontal velocity: dl=1,edim=2 means 1st 2D object
Inum=postint(fem,'u', 'dl',[1], 'edim',2)/fem.const.L;
% Show the result
disp(sprintf('Analytical flux = %g',-2/(fem.const.n+2)));
disp(sprintf('Numerical flux  = %g',Inum));