clc
clear
close all

% Space dimensions
fem.sdim={'x','y','z'};

% Geometry
g1=block3(1,1,1,'corner',zeros(1,3),[0 0 1],0);

clear s f c p
objs={g1};
name={'BLK1'};
s.objs=objs;
s.name=name;

objs={};
name={};
f.objs=objs;
f.name=name;

objs={};
name={};
c.objs=objs;
c.name=name;

objs={};
name={};
p.objs=objs;
p.name=name;

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

% Initialize mesh
fem.mesh=meshinit(fem,...
	'Out',    {'mesh'},...
	'jiggle', 'on',...
	'Hcurve', 0.40000000000000002,...
	'Hcutoff',0.01,...
	'Hgrad',  1.3999999999999999,...
	'Hmaxfact',1);

% Dimension
fem.dim={{'u'}};

% Boundary conditions
clear bnd
bnd.q={{{'0'}}};
bnd.g={{{'0'}}};
bnd.h={{{'0'}}};
bnd.r={{{'0'}}};
bnd.var={};
bnd.ind=ones(1,6);
fem.bnd=bnd;

% PDE coefficients
clear equ
equ.da={{{'1./sqrt(ux.^2+uy.^2+uz.^2)'}}};
equ.c={{{'1./sqrt(ux.^2+uy.^2+uz.^2)'}}};
equ.al={{{'0';'0';'0'}}};
equ.ga={{{'0';'0';'0'}}};
equ.be={{{'0';'0';'0'}}};
equ.a={{{'0'}}};
equ.f={{{'0'}}};
equ.var={'absux','sqrt(ux.^2+uy.^2+uz.^2)','abscu1x', ...
'sqrt(cu1x.^2+cu1y.^2+cu1z.^2)'};
equ.ind=1;
fem.equ=equ;

% Evaluate initial condition
fem.init=asseminit(fem,...
	'context','local',...
	'init',   struct('sd',{{{{'-((z-0.5).^2+(x-0.5).^2+(y-0.5).^2).^0.5+0.4'}}}}, ...
'ind',{1}));

% Solve dynamic problem
fem.sol=femtime(fem,...
	'tlist',  0:0.01:0.04,...
	'atol',   0.001,...
	'rtol',   0.01,...
	'jacobian','equ',...
	'mass',   'full',...
	'ode',    'ode15s',...
	'odeopt', struct('InitialStep',{[]},'MaxOrder',{5},'MaxStep',{[]}),...
	'out',    'sol',...
	'stop',   'on',...
	'report', 'on',...
	'context','local',...
	'sd',     'off',...
	'Epoint', 'gauss2',...
	'Nullfun','flnullorth',...
	'Tpoint', 'gauss2',...
	'Solcomp',1);

% Plot solution
figure(1)
postplot(fem,...
	'context','local',...
	'isodata',{'u','cont','on'},...
	'isofacestyle','none',...
	'isoedgestyle','interp',...
	'isolevels',0,...
	'isostyle','color',...
	'isomaxmin','off',...
	'isobar', 'on',...
	'isomap', 'bone',...
	'lightmodel','phong',...
	'lightreflection','default',...
	'transparency',1,...
	'view',   [-37.499999999999993 29.999999999999996],...
	'title',  'Time=0',...
	'solnum', 1,...
	'renderer','zbuffer',...
	'axisequal','on',...
	'axis',   [0 1 0 1 0 1],...
	'axisvisible','on',...
	'grid',   'off',...
	'campos', [-4.0657107175654046 -5.4501500521842638 4.8301270189221928],...
	'camprojection','orthographic',...
	'camtarget',0.5*ones(1,3),...
	'camup',  [0 0 1],...
	'camva',  14.437943202024892,...
	'axispos',[-0.024999999999999967 -0.05299999999999995 1.085 1.1409999999999998],...
	'scenelight','off',...
	'camlight','off',...
	'scenelightpos',[-4.0657107175654046 -5.4501500521842638 4.8301270189221928])

figure(2)
postplot(fem,...
	'context','local',...
	'isodata',{'u','cont','on'},...
	'isofacestyle','none',...
	'isoedgestyle','interp',...
	'isolevels',0,...
	'isostyle','color',...
	'isomaxmin','off',...
	'isobar', 'on',...
	'isomap', 'bone',...
	'lightmodel','phong',...
	'lightreflection','default',...
	'transparency',1,...
	'view',   [-37.499999999999993 29.999999999999996],...
	'title',  'Time=0.01',...
	'solnum', 2,...
	'renderer','zbuffer',...
	'axisequal','on',...
	'axis',   [0 1 0 1 0 1],...
	'axisvisible','on',...
	'grid',   'off',...
	'campos', [-4.0657107175654046 -5.4501500521842638 4.8301270189221928],...
	'camprojection','orthographic',...
	'camtarget',0.5*ones(1,3),...
	'camup',  [0 0 1],...
	'camva',  14.437943202024892,...
	'axispos',[-0.024999999999999967 -0.05299999999999995 1.085 1.1409999999999998],...
	'scenelight','off',...
	'camlight','off',...
	'scenelightpos',[-4.0657107175654046 -5.4501500521842638 4.8301270189221928])

figure(3)
postplot(fem,...
	'context','local',...
	'isodata',{'u','cont','on'},...
	'isofacestyle','none',...
	'isoedgestyle','interp',...
	'isolevels',0,...
	'isostyle','color',...
	'isomaxmin','off',...
	'isobar', 'on',...
	'isomap', 'bone',...
	'lightmodel','phong',...
	'lightreflection','default',...
	'transparency',1,...
	'view',   [-37.499999999999993 29.999999999999996],...
	'title',  'Time=0.02',...
	'solnum', 3,...
	'renderer','zbuffer',...
	'axisequal','on',...
	'axis',   [0 1 0 1 0 1],...
	'axisvisible','on',...
	'grid',   'off',...
	'campos', [-4.0657107175654046 -5.4501500521842638 4.8301270189221928],...
	'camprojection','orthographic',...
	'camtarget',0.5*ones(1,3),...
	'camup',  [0 0 1],...
	'camva',  14.437943202024892,...
	'axispos',[-0.024999999999999967 -0.05299999999999995 1.085 1.1409999999999998],...
	'scenelight','off',...
	'camlight','off',...
	'scenelightpos',[-4.0657107175654046 -5.4501500521842638 4.8301270189221928])

figure(3)
postplot(fem,...
	'context','local',...
	'isodata',{'u','cont','on'},...
	'isofacestyle','none',...
	'isoedgestyle','interp',...
	'isolevels',0,...
	'isostyle','color',...
	'isomaxmin','off',...
	'isobar', 'on',...
	'isomap', 'bone',...
	'lightmodel','phong',...
	'lightreflection','default',...
	'transparency',1,...
	'view',   [-37.499999999999993 29.999999999999996],...
	'title',  'Time=0.03',...
	'solnum', 4,...
	'renderer','zbuffer',...
	'axisequal','on',...
	'axis',   [0 1 0 1 0 1],...
	'axisvisible','on',...
	'grid',   'off',...
	'campos', [-4.0657107175654046 -5.4501500521842638 4.8301270189221928],...
	'camprojection','orthographic',...
	'camtarget',0.5*ones(1,3),...
	'camup',  [0 0 1],...
	'camva',  14.437943202024892,...
	'axispos',[-0.024999999999999967 -0.05299999999999995 1.085 1.1409999999999998],...
	'scenelight','off',...
	'camlight','off',...
	'scenelightpos',[-4.0657107175654046 -5.4501500521842638 4.8301270189221928])

figure(4)
postplot(fem,...
	'context','local',...
	'isodata',{'u','cont','on'},...
	'isofacestyle','none',...
	'isoedgestyle','interp',...
	'isolevels',0,...
	'isostyle','color',...
	'isomaxmin','off',...
	'isobar', 'on',...
	'isomap', 'bone',...
	'lightmodel','phong',...
	'lightreflection','default',...
	'transparency',1,...
	'view',   [-37.499999999999993 29.999999999999996],...
	'title',  'Time=0.04',...
	'solnum', 5,...
	'renderer','zbuffer',...
	'axisequal','on',...
	'axis',   [0 1 0 1 0 1],...
	'axisvisible','on',...
	'grid',   'off',...
	'campos', [-4.0657107175654046 -5.4501500521842638 4.8301270189221928],...
	'camprojection','orthographic',...
	'camtarget',0.5*ones(1,3),...
	'camup',  [0 0 1],...
	'camva',  14.437943202024892,...
	'axispos',[-0.024999999999999967 -0.05299999999999995 1.085 1.1409999999999998],...
	'scenelight','off',...
	'camlight','off',...
	'scenelightpos',[-4.0657107175654046 -5.4501500521842638 4.8301270189221928])