Integrating gradient of harmonic function over 2D domain from PDE Solver

I am trying to evaluate over a specific domain (described below), where is a harmonic function solved using Matlab's PDE solver.

• In particular, how can I access the values of from the PDE solver app?
• Is there a method to evalute automatically with Matlab?
• Or should I approach this problem by summing: where is the area of the triangles generated by the mesh function from PDE solver.

Thanks in advance.

% This script is written and read by pdetool and should NOT be edited.
% There are two recommended alternatives:
% 1) Export the required variables from pdetool and create a MATLAB script
%    to perform operations on these.
% 2) Define the problem completely using a MATLAB script. See
%    http://www.mathworks.com/help/pde/examples/index.html for examples
%    of this approach.



function pdemodel
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',1);
set(ax,'DataAspectRatio',[1 1 1]);
set(ax,'PlotBoxAspectRatio',[665.60000000000002 450 96.428571428571445]);
set(ax,'XLimMode','auto');
set(ax,'YLimMode','auto');
set(ax,'XTickMode','auto');
set(ax,'YTickMode','auto');
pdetool('gridon','on');

% Geometry description:
% e = 0.005, 0.5, 0.125
% L = 0.5, 1, 2, 4, 8 


e = 0.001;
L = 20;

pderect([-L/2 L/2 1 e],'R1');
pderect([-L/2 L/2 -1 -e],'R2');
pdecirc(L/2,0,1,'C1');
pdecirc(-L/2,0,1,'C2');
pdecirc(L/2,0,e,'C3');
pdecirc(-L/2,0,e,'C4');
pderect([-L/2 L/2 e -e],'R3');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','(R1+R2+C1+C2) - (C3 + R3 + C4)')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(24,...
'dir',...
1,...
'1',...
'0')
pdesetbd(23,...
'dir',...
1,...
'1',...
'0')
pdesetbd(22,...
'dir',...
1,...
'1',...
'0')
pdesetbd(21,...
'dir',...
1,...
'1',...
'0')
pdesetbd(20,...
'dir',...
1,...
'1',...
'1')
pdesetbd(17,...
'dir',...
1,...
'1',...
'1')
pdesetbd(15,...
'dir',...
1,...
'1',...
'1')
pdesetbd(14,...
'dir',...
1,...
'1',...
'1')
pdesetbd(12,...
'dir',...
1,...
'1',...
'0')
pdesetbd(11,...
'dir',...
1,...
'1',...
'0')
pdesetbd(10,...
'dir',...
1,...
'1',...
'0')
pdesetbd(9,...
'dir',...
1,...
'1',...
'0')
pdesetbd(8,...
'dir',...
1,...
'1',...
'0')
pdesetbd(7,...
'dir',...
1,...
'1',...
'0')
pdesetbd(2,...
'dir',...
1,...
'1',...
'1')
pdesetbd(1,...
'dir',...
1,...
'1',...
'1')

% Mesh generation:
setappdata(pde_fig,'Hgrad',1.3);
setappdata(pde_fig,'refinemethod','regular');
setappdata(pde_fig,'jiggle',char('on','mean',''));
setappdata(pde_fig,'MesherVersion','preR2013a');
pdetool('initmesh')
pdetool('refine')

% PDE coefficients:
pdeseteq(1,...
'1.0',...
'0.0',...
'0',...
'1.0',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setappdata(pde_fig,'currparam',...
['1.0';...
'0.0';...
'0  ';...
'1.0'])

% Solve parameters:
setappdata(pde_fig,'solveparam',...
char('0','3216','10','pdeadworst',...
'0.5','longest','0','1E-4','','fixed','Inf'))

% Plotflags and user data strings:
setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 1 0 0 1]);
setappdata(pde_fig,'colstring','');
setappdata(pde_fig,'arrowstring','');
setappdata(pde_fig,'deformstring','');
setappdata(pde_fig,'heightstring','');

% Solve PDE:
pdetool('solve');