% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5a (COMSOL 3.5.0.606, $Date: 2009/04/29 09:11:29 $)
% 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.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 606;
vrsn.rcs = '$Name: v35ap $';
vrsn.date = '$Date: 2009/04/29 09:11:29 $';
fem.version = vrsn;


% Geometry
g1=sphere3('1e-4','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g2=sphere3('2e-4','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g3=geomcomp({g1,g2},'ns',{'g1','g2'},'sf','g1+g2','face','none','edge','all');
g4=block3('1e-3','1e-3','1e-3','base','center','pos',{'0','0','-.5e-4'},'axis',{'0','0','1'},'rot','0');
g4=move(g4,[0,0,-.5e-4]);
g5=block3('1e-3','1e-3','1e-3','base','center','pos',{'0','0','-.5e-3'},'axis',{'0','0','1'},'rot','0');
g6=geomcomp({g3,g5},'ns',{'g3','g5'},'sf','g3*g5','face','none','edge','all');
parr={point3(0,0,0)};
g7=geomcoerce('point',parr);

% Analyzed geometry
clear p s
p.objs={g7};
p.name={'PT1'};
p.tags={'g7'};

s.objs={g6};
s.name={'CO2'};
s.tags={'g6'};

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

% flbinaryfile='nickel_model.mphm';

% Constants
fem.const = {'r0','4e-7', ...
  't0','1e-9', ...
  'r1','1e-4', ...
  'r2','2e-4', ...
  'I0','1e12'};

% Geometry
% fem.geom = flbinary('g1','geom',flbinaryfile);

% Application mode 1
clear appl
appl.mode.class = 'SmeSolid3';
appl.module = 'SME';
appl.gporder = 4;
appl.cporder = 2;
appl.assignsuffix = '_smsld';
clear prop
prop.analysis='time';
appl.prop = prop;
clear equ
equ.nu = 0.35;
equ.Tflag = 1;
equ.rho = 8912;
equ.betadK = {'((r-r1)/(r2-r1))^3*1e-6',0};
equ.dampingtype = {'Rayleigh','nodamping'};
 equ.D = {{1.6360e11,0.9230e11,0.6792e11,0,0,0;0.9230e11, ...
   1.6360e11,0.6792e11,0,0,0;0.6792e11,0.6792e11,1.8520e11, ...
   0,0,0;0,0,0,0.3565e11,0,0;0,0,0,0,0.4705e11, ...
   0;0,0,0,0,0,0.4705e11}};
equ.D = C_cell;
equ.E = 110e9;
equ.Tempref = 273.15;
equ.materialmodel = 'aniso';
equ.Temp = 'T';
equ.alpha = 17e-6;
equ.name = {'','default'};
equ.alphavector = {{1.34e-5;1.34e-5;1.34e-5}};
equ.ind = [1,2];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.assignsuffix = '_ht';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm4'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.q0 = {0,0,'q'};
bnd.type = {'q0','cont','q'};
bnd.ind = [1,1,1,2,3,2,1,2,1,2];
appl.bnd = bnd;
clear equ
equ.C = '445[J/(kg*K)]';
equ.k = '90.7[W/(m*K)]';
equ.name = 'default';
equ.rho = '8912[kg/m^3]';
equ.ind = [1,1];
appl.equ = equ;
fem.appl{2} = appl;
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

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

% Subdomain expressions
equ.expr = {'r',{'sqrt(x^2+y^2+z^2)',''}};
fem.equ = equ;

% Boundary settings
clear bnd
bnd.ind = [1,1,1,1,2,1,1,1,1,1];

% Boundary expressions
bnd.expr = {'rxy',{'','sqrt(x^2+y^2)'}, ...
  'q',{'','I0*t*exp(-t/t0)*exp(-rxy^2/r0^2)/t0'}};
fem.bnd = bnd;

% Library materials
 clear lib
 lib.mat{1}.name='Ni';
 lib.mat{1}.varname='mat1';
 lib.mat{1}.variables.nu='0.31';
 lib.mat{1}.variables.E='219e9[Pa]';
 lib.mat{1}.variables.sigma='13.8e6[S/m]';
 lib.mat{1}.variables.alpha='13.4e-6[1/K]';
 lib.mat{1}.variables.C='445[J/(kg*K)]';
 lib.mat{1}.variables.rho='8900[kg/m^3]';
 lib.mat{1}.variables.k='90.7[W/(m*K)]';
 
 fem.lib = lib;

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

% Multiphysics
fem=multiphysics(fem);

% Initialise the mesh 
fem.mesh=meshinit(fem, ...
    'hauto',mesh_level, ...
    'hmaxvtx',[7,1e-7], ...
    'hgradvtx',[7,10], ...
    'hmaxfac',[5,4.e-6], ...
    'hgradfac',[5,10]);

fem.xmesh = meshextend(fem);

% Evaluate initial value
 init = asseminit(fem,'blocksize','auto');

% Update model
 fem.sol = init;

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

% Solve the model
fem.sol=femtime(fem, ...
    'solcomp',{'w','T','v','u'}, ...
    'outcomp',{'w','v','T','u'}, ...
    'blocksize','auto', ...
    'odesolver','genalpha', ...
    'tlist',times, ...
    'tout','tlist', ...
    'atol',{'1e-16'}, ...
    'incrdelay','off', ...
    'nlsolver','manual', ...
    'linsolver','spooles', ...
    'ntolfact',1, ...
    'maxiter',4, ...
    'dtech','const', ...
    'damp',1.0, ...
    'jtech','minimal');

% Plot solution
 postplot(fem, ...
          'tridata',{'w','cont','internal','unit','m'}, ...
          'tridlim',[-5e-12 5e-12], ...
          'trimap','Rainbow', ...
          'title','Time=0    Boundary: z-displacement [m]', ...
          'geom','off', ...
          'grid','on', ...
          'campos',[-0.0021924694954226307,9.008658625091105E-4,0.001798428853438161], ...
          'camtarget',[-3.2522338210578045E-5,-8.03394678989124E-6,-7.467791340396796E-5], ...
          'camup',[0.6096405526049752,-0.1537799791553065,0.7776182319303567], ...
          'camva',14.673518253685748);