% FEMLAB Model M-file
% Generated 10-May-2003 19:36:55 by FEMLAB 2.3.0.148.

flclear fem
% FEMLAB Version
clear vrsn;
vrsn.name='FEMLAB 2.3';
vrsn.major=0;
vrsn.build=148;
fem.version=vrsn;

% Recorded command sequence

% New geometry 1
fem.sdim={'x','y'};

% Geometry
clear s c p
R1=rect2(0,1,0,0.5,0);
objs={R1};
names={'R1'};
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);

clear appl

% Application mode 1
appl{1}.mode=flheat2d('dim',{'u','u_t'},'sdim',{'x','y'},'submode','std', ...
'tdiff','on');
appl{1}.dim={'u','u_t'};
appl{1}.form='coefficient';
appl{1}.border='off';
appl{1}.name='c1';
appl{1}.var={};
appl{1}.assign={'abscu1x';'abscu1x';'absga1x';'absga1x';'absux';'absux'};
appl{1}.elemdefault='Lag2';
appl{1}.shape={'shlag(2,''u'')'};
appl{1}.sshape=2;
appl{1}.equ.da={{{'1'}}};
appl{1}.equ.c={{{'1'}}};
appl{1}.equ.f={{{'1'}}};
appl{1}.equ.gporder={{4}};
appl{1}.equ.cporder={{2}};
appl{1}.equ.shape={1};
appl{1}.equ.init={{{'0'}}};
appl{1}.equ.usage={1};
appl{1}.equ.ind=1;
appl{1}.bnd.q={{{'0'}}};
appl{1}.bnd.g={{{'0'}}};
appl{1}.bnd.h={{{'1'}}};
appl{1}.bnd.r={{{'0'}}};
appl{1}.bnd.type={'dir'};
appl{1}.bnd.gporder={{0}};
appl{1}.bnd.cporder={{0}};
appl{1}.bnd.shape={0};
appl{1}.bnd.ind=ones(1,4);

fem.appl=appl;

% Initialize mesh
fem.mesh=meshinit(fem,...
	'Out',    {'mesh'},...
	'jiggle', 'off',...
	'Hcurve', 0.29999999999999999,...
	'Hgrad',  1.3,...
	'Hmax',   {0.5,zeros(1,0),zeros(1,0),zeros(1,0)},...
	'Hnum',   {[],zeros(1,0)},...
	'Hpnt',   {10,zeros(1,0)});

% Geometry
clear s c p
R1=rect2(-1,1,0,1,0);
objs={R1};
names={'R1'};
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', 'off',...
	'Hcurve', 0.29999999999999999,...
	'Hgrad',  1.3,...
	'Hmax',   {0.5,zeros(1,0),zeros(1,0),zeros(1,0)},...
	'Hnum',   {[],zeros(1,0)},...
	'Hpnt',   {10,zeros(1,0)});

% Refine mesh
fem.mesh=meshrefine(fem,...
	'out',    {'mesh'},...
	'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem,...
	'out',    {'mesh'},...
	'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem,...
	'out',    {'mesh'},...
	'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem,...
	'out',    {'mesh'},...
	'rmethod','regular');

% Differentiation rules
fem.rules={};

% Problem form
fem.outform='coefficient';

% Differentiation
fem.diff={'expr'};

% Differentiation simplification
fem.simplify='on';

% Boundary conditions
clear bnd
bnd.q={{{'0'}},{{'0'}},{{'0'}}};
bnd.g={{{'0'}},{{'0'}},{{'0'}}};
bnd.h={{{'1'}},{{'1'}},{{'1'}}};
bnd.r={{{'5'}},{{'0'}},{{'0'}}};
bnd.type={'dir','neu','dir'};
bnd.gporder={{0},{0},{0}};
bnd.cporder={{0},{0},{0}};
bnd.shape={0,0,0};
bnd.ind=[1 2 3 3];
fem.appl{1}.bnd=bnd;

% PDE coefficients
clear equ
equ.da={{{'1'}}};
equ.c={{{'1'}}};
equ.f={{{'0'}}};
equ.gporder={{4}};
equ.cporder={{2}};
equ.shape={1};
equ.init={{{'0'}}};
equ.usage={1};
equ.ind=1;
fem.appl{1}.equ=equ;

% Internal borders
fem.appl{1}.border='off';

% Shape functions
fem.appl{1}.shape={'shlag(2,''u'')'};

% Geometry element order
fem.appl{1}.sshape=2;

% Define constants
fem.const={};

% Multiphysics
fem=multiphysics(fem);

% Extend the mesh
fem.xmesh=meshextend(fem,'context','local','cplbndeq','on','cplbndsh','on');

% Evaluate initial condition
init=asseminit(fem,...
	'context','local',...
	'init',   fem.xmesh.eleminit);

% Solve dynamic problem
fem.sol=femtime(fem,...
	'tlist',  0:0.001:1,...
	'atol',   1e-006,...
	'rtol',   1e-006,...
	'jacobian','equ',...
	'mass',   'full',...
	'ode',    'ode15s',...
	'odeopt', struct('InitialStep',{[]},'MaxOrder',{5},'MaxStep',{[]}),...
	'out',    'sol',...
	'stop',   'on',...
	'init',   init,...
	'report', 'on',...
	'timeind','auto',...
	'context','local',...
	'sd',     'off',...
	'nullfun','flnullorth',...
	'blocksize',5000,...
	'solcomp',{'u'},...
	'linsolver','matlab',...
	'uscale', 'auto');

% Plot solution
postplot(fem,...
	'geomnum',1,...
	'context','local',...
	'tridata',{'u','cont','internal'},...
	'trifacestyle','interp',...
	'triedgestyle','none',...
	'trimap', 'jet',...
	'trimaxmin','off',...
	'tribar', 'off',...
	'geom',   'on',...
	'geomcol','bginv',...
	'refine', 3,...
	'contorder',2,...
	'phase',  0,...
    'solnum', 425,...
	'title',  '',...
	'renderer','zbuffer',...
    'contdata','u',...
    'contmap','cool',...
    'contlevels',0:0.2:5,...
	'axisvisible','on')
set(gca,'FontName','Times','FontSize',16,'Box','on')
xlabel('{\itx}','FontSize',21)
ylabel('{\ity}','FontSize',21)
axis equal
axis tight
print -deps2c u424color.eps

uu=postinterp(fem,'u',[0;0],'solnum',1:length(fem.sol.tlist));

figure(2);
plot(fem.sol.tlist,uu,'k-','LineWidth',2);
set(gca,'FontName','Times','FontSize',21,'Box','on')
set(gca,'XTickMode','manual','XTick',[0 0.424 1],...
        'YTickMode','manual','YTick',[0 1 1.25]);
ylabel('{\itu}({\itt},0,0)','FontWeight','normal','FontSize',28);
xlabel('time{} {\itt}','FontWeight','normal','FontSize',28);
axis([-0.008 1.008 -0.01 1.26]);
grid on
print -deps2c u00.eps

t1 = interp1(uu,fem.sol.tlist,1,'spline')