function fem = getfem (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=rect2(1,1,'base','corner','pos',[0,0]);

% Analyzed geometry
clear s
s.objs={g1};
s.name={'R1'};
s.tags={'g1'};

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

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',9);

% 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.assignsuffix = '_c';
clear prop
prop.elemdefault='Lag1';
appl.prop = prop;
clear bnd
bnd.type = 'dir';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
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',[-0.2519372693726938,1.2519372693726938,-0.05,1.05,-1,1]);

% Plot solution
% postplot(fem, ...
%          'tridata',{'u','cont','internal'}, ...
%          'triz','u', ...
%          'trimap','jet(1024)', ...
%          'title','Surface: u   Height: u', ...
%          'grid','on');