% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5a (COMSOL 3.5.0.603, $Date: 2008/12/03 17:02:19 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 603;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/12/03 17:02:19 $';
fem.version = vrsn;

% Constants
fem.const = {'Pin','2989.0 [N/m^2]', ...
  'Pout','0 [N/m^2]', ...
  'rho_water','995.7 [kg/m^3]', ...
  'eta_water','0.798 [Pa *s]', ...
  'm_tube','1.5e6 [Pa]', ...
  'v_tube','20*m_tube'};

% Geometry
g1=cylinder3('.003175','.0762','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g2=cylinder3('.0028702','.0762','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g3=geomcomp({g1},'ns',{'g1'},'sf','g1','face','none','edge','all');
g4=geomcomp({g3},'ns',{'g3'},'sf','g3','face','none','edge','all');
g5=geomcomp({g4,g2},'ns',{'g4','g2'},'sf','g4-g2','face','none','edge','all');

% Analyzed geometry
clear s
s.objs={g5};
s.name={'CO1'};
s.tags={'g5'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);
g6=cylinder3('.0028702','.0762','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');

% Analyzed geometry
clear s
s.objs={g5,g6};
s.name={'CO1','CYL1'};
s.tags={'g5','g6'};

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

% Constants
fem.const = {'Pin','2989.0 [N/m^2]', ...
  'Pout','0 [N/m^2]', ...
  'rho_water','995.7 [kg/m^3]', ...
  'eta_water','0.798 [Pa *s]', ...
  'm_tube','1.5e6 [Pa]', ...
  'v_tube','20*m_tube', ...
  'rho_tube','920 [kg/m^3]'};

% Constants
fem.const = {'Pin','2989.0 [N/m^2]', ...
  'Pout','2000 [N/m^2]', ...
  'rho_water','995.7 [kg/m^3]', ...
  'eta_water','0.798 [Pa *s]', ...
  'm_tube','1.5e6 [Pa]', ...
  'v_tube','20*m_tube', ...
  'rho_tube','920 [kg/m^3]'};

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4, ...
                  'hmaxsub',[2,1e-3]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'FlNavierStokes';
appl.gporder = 2;
appl.cporder = 1;
appl.assignsuffix = '_ns';
clear prop
prop.elemdefault='Lag1';
appl.prop = prop;
clear bnd
bnd.p0 = {0,0,'Pout*amp','Pin'};
bnd.type = {'int','walltype','outlet','inlet'};
bnd.outtype = {'p','p','p0','p'};
bnd.intype = {'uv','uv','uv','p'};
bnd.ind = [1,1,1,1,2,2,3,4,1,2,2,1];
appl.bnd = bnd;
clear equ
equ.eta = {1,'eta_water'};
equ.rho = {1,'rho_water'};
equ.usage = {0,1};
equ.ind = [1,2];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'SmeSolid3';
appl.dim = {'u2','v2','w2','p2'};
appl.module = 'SME';
appl.shape = {'shlag(1,''u2'')','shlag(1,''v2'')','shlag(1,''w2'')','shlag(2,''u2'')','shlag(2,''v2'')','shlag(2,''w2'')'};
appl.gporder = {2,4};
appl.cporder = {1,2};
appl.assignsuffix = '_smsld';
clear prop
prop.analysis='quasi';
prop.largedef='on';
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm5','lm6','lm7'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.name = {'Outer','Ends','Inner'};
bnd.Fz = {0,0,'-T_z_ns'};
bnd.Fx = {0,0,'-T_x_ns'};
bnd.Fy = {0,0,'-T_y_ns'};
bnd.constrcond = {'free','roller','free'};
bnd.ind = [1,1,2,2,3,3,1,1,1,3,3,1];
appl.bnd = bnd;
clear equ
equ.kappa = {'v_tube',1e10};
equ.gporder = {1,2};
equ.rho = {'rho_tube',7850};
equ.cporder = {1,2};
equ.mu = {'m_tube',8e5};
equ.materialmodel = {'hyper','iso'};
equ.shape = {[1;2;3],[4;5;6]};
equ.alpha = {1.5e-4,1.2e-5};
equ.usage = {1,0};
equ.ind = [1,2];
appl.equ = equ;
fem.appl{2} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'amp','(t<=0.5)*sin(pi*t[1/s])+ (t>0.5)*(1.5-0.5*cos(-2*pi*(0.5-t[1/s])))'};

% Descriptions
clear descr
descr.const= {'Pin','inlet pressure','rho_water','density of water','m_tube','Neo-hookean hyperelastic behavior, mu coefficient','v_tube','bulk modulus of tube','Pout','outlet pressure','rho_tube','density of tube','eta_water','dynamic viscosity of water'};
fem.descr = descr;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

% Multiphysics
fem=multiphysics(fem);

% Generate GMG mesh cases
fem=meshcaseadd(fem,'mgauto','anyshape');

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

% Solve problem
fem.sol=femtime(fem, ...
                'symmetric','off', ...
                'solcomp',{'w','v','u','p'}, ...
                'outcomp',{'w','v','u','u2','p','v2','w2'}, ...
                'blocksize','auto', ...
                'odesolver','genalpha', ...
                'tlist',[colon(0,0.15,1)], ...
                'estrat',1, ...
                'tout','tlist', ...
                'tsteps','manual', ...
                'rhoinf',0.5, ...
                'timestep',.01, ...
                'nlsolver','manual', ...
                'linsolver','bicgstab', ...
                'rhob',20, ...
                'maxlinit',100, ...
                'prefuntype','right', ...
                'prefun','gmg', ...
                'prepar',{'presmooth','vanka','presmoothpar',{'iter',1,'vankasolv','direct','seconditer',3,'relax',0.8},'postsmooth','vanka','postsmoothpar',{'iter',1,'vankasolv','direct','seconditer',3,'relax',0.8},'csolver','pardiso'}, ...
                'ntolfact',1, ...
                'maxiter',5, ...
                'dtech','const', ...
                'damp',1.0, ...
                'jtech','once', ...
                'uscale','none', ...
                'mcase',[0 1]);

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

% Plot solution
postplot(fem, ...
         'slicedata',{'U_ns','cont','internal','unit','m/s'}, ...
         'slicexspacing',5, ...
         'sliceyspacing',0, ...
         'slicezspacing',0, ...
         'slicemap','Rainbow', ...
         'solnum','end', ...
         'title','Time=0.9    Slice: Velocity field [m/s]', ...
         'grid','on', ...
         'campos',[0.05784795862965936,0.23499522510413187,-0.2595709359673219], ...
         'camtarget',[0,0,0.038100000470876694], ...
         'camup',[0.8852145939882246,-0.433064763329723,-0.16985297569496324], ...
         'camva',11.793824864599486);

% Remove generated GMG mesh cases
fem=meshcasedel(fem,[1]);

return; %%%%%%%%%%% END OF HISTORY %%%%%%%%%%%%%%%%%%%%%%%
%                                                        %
% M-file below is used when resetting the model.         %
%                                                        %
%%%%%%%%%%%%%% START OF RESET HISTORY %%%%%%%%%%%%%%%%%%%%

% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5a (COMSOL 3.5.0.603, $Date: 2008/12/03 17:02:19 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 603;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/12/03 17:02:19 $';
fem.version = vrsn;

% Constants
fem.const = {'Pin','2989.0 [N/m^2]', ...
  'Pout','0 [N/m^2]', ...
  'rho_water','995.7 [kg/m^3]', ...
  'eta_water','0.798 [Pa *s]', ...
  'm_tube','1.5e6 [Pa]', ...
  'v_tube','20*m_tube'};

% Geometry
g1=cylinder3('.003175','.0762','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g2=cylinder3('.0028702','.0762','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g3=geomcomp({g1},'ns',{'g1'},'sf','g1','face','none','edge','all');
g4=geomcomp({g3},'ns',{'g3'},'sf','g3','face','none','edge','all');
g5=geomcomp({g4,g2},'ns',{'g4','g2'},'sf','g4-g2','face','none','edge','all');

% Analyzed geometry
clear s
s.objs={g5};
s.name={'CO1'};
s.tags={'g5'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);
g6=cylinder3('.0028702','.0762','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');

% Analyzed geometry
clear s
s.objs={g5,g6};
s.name={'CO1','CYL1'};
s.tags={'g5','g6'};

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

% Constants
fem.const = {'Pin','2989.0 [N/m^2]', ...
  'Pout','2000 [N/m^2]', ...
  'rho_water','995.7 [kg/m^3]', ...
  'eta_water','0.798 [Pa *s]', ...
  'm_tube','1.5e6 [Pa]', ...
  'v_tube','20*m_tube', ...
  'rho_tube','920 [kg/m^3]'};