nrefines = 1;

% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5a (COMSOL 3.5.0.603, $Date: 2008/12/03 17:02:19 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 603;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/12/03 17:02:19 $';
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 i = 1 : nrefines
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 = '(-2*pi^2)*(cos(2*pi*x)*sin(pi*y)^2+sin(pi*x)^2*cos(2*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;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

% Multiphysics
fem=multiphysics(fem);

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

% Solve problem
tic
fem.sol=femstatic(fem, ...
                  'solcomp',{'u'}, ...
                  'outcomp',{'u'}, ...
                  'blocksize','auto', ...
                  'linsolver','pardiso');
toc
% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'u','cont','internal'}, ...
         'trimap','Rainbow', ...
         'title','Surface: u', ...
         'axis',[-0.3835106382978725,1.3835106382978726,-0.05,1.05]);

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

print -djpeg100 poisson2d.jpg

% Integrate
I1=postint(fem,'(u-sin(pi*x)^2*sin(pi*y)^2)^2', ...
           'unit','', ...
           'recover','off', ...
           'dl',1);

errnorm = sqrt(I1);