function fem = getfem_cos (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(2,2,'base','center','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, ...
                  '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.assignsuffix = '_c';
clear prop
prop.elemdefault='Lag1';
appl.prop = prop;
clear bnd
bnd.type = 'dir';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.f = '(pi^2/2)*(cos(pi*x)*(cos(pi*y/2))^2+(cos(pi*x/2))^2*cos(pi*y))';
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'});