function fem = getfem_del_disk (nref)

% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.3 (COMSOL 3.3.0.405, $Date: 2006/08/31 18:03:47 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.3';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 405;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2006/08/31 18:03:47 $';
fem.version = vrsn;

% Geometry
g1=ellip2(1,1,'base','center','pos',[0,0]);
parr={point2(0,0)};
g2=geomcoerce('point',parr);

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

s.objs={g1};
s.name={'E1'};
s.tags={'g1'};

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

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hmax',[2], ...
                  'hmaxfact',5, ...
                  'hgrad',2, ...
                  'hcurve',1, ...
                  'hcutoff',0.05);

% Refine mesh
for nr = 1 : nref
  fem.mesh=meshrefine(fem, ...
                      'mcase',0, ...
                      'rmethod','regular');
end;

% (Default values are not included)

% Application mode 1
% clear appl
% appl.mode.class = 'FlPDEC';
% appl.shape = {'shlag(1,''u'')'};
% appl.gporder = {2,4};
% appl.cporder = 1;
% appl.assignsuffix = '_c';
% alternative way to implement element choice, most suitable to Lag2, etc.:
clear appl
appl.mode.class = 'FlPDEC';
appl.assignsuffix = '_c';
clear prop
prop.elemdefault='Lag1';
appl.prop = prop;
clear pnt
pnt.weak = {0,'u_test'};
pnt.ind = [1,1,2,1,1];
appl.pnt = pnt;
clear bnd
bnd.type = 'dir';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.f = 0;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Multiphysics
fem=multiphysics(fem);

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

% Solve problem
fem.sol=femstatic(fem, ...
                  'solcomp',{'u'}, ...
                  'outcomp',{'u'});

%% Save current fem structure for restart purposes
%fem0=fem;
%
%% Plot solution
%postplot(fem, ...
%         'tridata',{'u','cont','internal'}, ...
%         'trimap','jet(1024)', ...
%         'title','Surface: u', ...
%         'axis',[-1.453136531365314,1.453136531365314,-1.1,1.1,-1,1]);
%
%% Plot solution
%postplot(fem, ...
%         'tridata',{'u','cont','internal'}, ...
%         'triz','u', ...
%         'trimap','jet(1024)', ...
%         'title','Surface: u   Height: u', ...
%         'grid','on');
%
%% Integrate
%I1=postint(fem,'(u+log(sqrt(x^2+y^2))/(2*pi))^2', ...
%           'unit','', ...
%           'dl',[1]);