close all;
clear;

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

% Geometry
clear s c p
R1=rect2(0,2,0,2,0);
R2=rect2(0,2,8,10,0);
R3=rect2(8,10,8,10,0);
R4=rect2(8,10,0,2,0);
R5=rect2(2,8,2,8,0);
R6=rect2(0,2,2,8,0);
R7=rect2(2,8,8,10,0);
R8=rect2(8,10,2,8,0);
R9=rect2(2,8,0,2,0);
objs={R1,R2,R3,R4,R5,R6,R7,R8,R9};
names={'R1','R2','R3','R4','R5','R6','R7','R8','R9'};
s.objs=objs;
s.name=names;

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

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

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

% Initialize mesh
fem.mesh=meshinit(fem,...
	'Out',    {'mesh'},...
	'jiggle', 'mean',...
	'Hcurve', 0.29999999999999999,...
	'Hgrad',  1.3,...
	'Hpnt',   {10,[]});

fem.mesh = meshrefine(fem);

% 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,24);
fem.bnd=bnd;

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

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

% Solve dynamic problem
fem.sol=femtime(fem,...
	'tlist',  0:0.1:6,...
	'atol',   1e-8,...
	'rtol',   1e-6,...
	'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',...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
    'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=0 ',...
	'renderer','zbuffer',...
	'solnum', 1,...
	'axisvisible','on')
axis equal;

% Plot solution
figure(2);
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=0.5' ,...
	'renderer','zbuffer',...
	'solnum', 6,...
	'axisvisible','on')
axis equal;


% Plot solution
figure(3);
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=1',...
	'renderer','zbuffer',...
	'solnum', 11,...
	'axisvisible','on')
axis equal;


% Plot solution
figure(4)
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=1.1  ',...
	'renderer','zbuffer',...
	'solnum', 12,...
	'axisvisible','on')
axis equal;


% Plot solution
figure(5)
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=1.2',...
	'renderer','zbuffer',...
	'solnum', 13,...
	'axisvisible','on')
axis equal;


% Plot solution
figure(6)
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=1.3',...
	'renderer','zbuffer',...
	'solnum', 14,...
	'axisvisible','on')
axis equal;


% Plot solution
figure(7)
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=1.4',...
	'renderer','zbuffer',...
	'solnum', 15,...
	'axisvisible','on')
axis equal;


% Plot solution
figure(8)
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=2',...
	'renderer','zbuffer',...
	'solnum', 21,...
	'axisvisible','on')
axis equal;


% Plot solution
figure(9)
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=4',...
	'renderer','zbuffer',...
	'solnum', 41,...
	'axisvisible','on')
axis equal;

figure(10)
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=6',...
	'renderer','zbuffer',...
	'solnum', 61,...
	'axisvisible','on')
axis equal;