Computing electrostatic force in pde toolbox

I want to compute the electrostatic force on a conductor after solving the Poisson equation in the domain. The code to solve the Poisson equation is given below (This is based on the code given in a pde example: Finite Element Analysis of Electrostatically Actuated MEMS Device)

h = 1;
g0 = h;
b = 5*h;
bg = b;
LX =5*bg;
LY = 2*g0+h;
hmax = 0.2;
 x1  = [0;(LX-bg)/2;(LX-b)/2;LX/2;(LX+b)/2;(LX+bg)/2;LX]; % X-coord of init mesh
 y1  = [0;g0;g0+h/2;g0+h;LY];                             % Y-coord of init mesh

rect_domain = [2,6,x1(1), x1(2), x1(6) x1(7),  x1(7), x1(1),y1(1),y1(1),y1(1),y1(1),y1(5),y1(5)]'; % electrostatic domain


rect_movable = [3,4,x1(3)  x1(5), x1(5), x1(3),y1(2), y1(2),y1(4),y1(4)]';% Conductor region

rect_movable = [rect_movable;zeros(length(rect_domain)-length(rect_movable),1)]; % padding zeros to match with the size of rect_domain

gd = [rect_domain,rect_movable];

ns = char('rect_domain','rect_movable');
ns = ns';

sf = 'rect_domain-(rect_movable)';
dl = decsg(gd,sf,ns);

model = createpde;
geometryFromEdges(model,dl);

pdegplot(model,'EdgeLabels','on','FaceLabels','on')
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([0,LX,0, LY])
axis tight

V0 = 0;
V1 = 10;
applyBoundaryCondition(model,'dirichlet','Edge',[6],'u',V0); % ground electrode 
applyBoundaryCondition(model,'dirichlet','Edge',[3,4,8,9],'u',V1); % Conductor

specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',0);

% Generate mesh
generateMesh(model,'Hmax',hmax);
pdeplot(model)
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([0,LX,0, LY])
axis tight

% Poisson solver
results = solvepde(model);

u = results.NodalSolution;
figure
pdeplot(model,'XYData',results.NodalSolution,...
              'ColorMap','jet');

title('Electric Potential');
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([0,LX,0, LY])
           

In the example cited( Finite Element Analysis of Electrostatically Actuated MEMS Device ), They have used a function handle to compute the electrostatic pressure as shown below;

pressureFcn = @(location,state) - ...
                CalculateElectrostaticPressure(results,[],location);

The function used in the function handle has been given as;

function ePressure = ...
    CalculateElectrostaticPressure(elecResults,structResults,location)
% Function to compute electrostatic pressure.
% structuralBoundaryLoad is used to specify the pressure load
% on the movable electrode.
% Inputs:
% elecResults: Electrostatic FEA results
% structResults: Mechanical FEA results (optional)
% location: The x,y coordinate where pressure is obtained
%
% Output:
% ePressure: Electrostatic pressure at location
%
% location.x : The x-coordinate of the points
% location.y : The y-coordinate of the points
xq = location.x;
yq = location.y;
% Compute the magnitude of the electric field from the potential
% difference.
[gradx,grady] = evaluateGradient(elecResults,xq,yq);
absE = sqrt(gradx.^2 + grady.^2);
% The permittivity of vacuum is 8.854*10^-12 farad/meter.
epsilon0 = 8.854e-12;

% Compute the magnitude of the electric flux density.
absD0 = epsilon0*absE;
absD = absD0;

% If structResults (deformation) is available,
% update the charge density along the movable electrode.
if ~isempty(structResults)
    % Displacement of the movable electrode at position x
    intrpDisp = interpolateDisplacement(structResults,xq,yq);
    vdisp = abs(intrpDisp.uy);
    G = 2e-6; % Gap 2 micron
    absD = absD0.*G./(G-vdisp);
end

% Compute the electrostatic pressure.
ePressure = absD.^2/(2*epsilon0);

end

Here, I have two queries;

• What are the arguments (location,state) passed into the function CalculateElectrostaticPressure. I have to compute the electrostatic pressure on the surface of the conductor) - I have little or no clue about how this function handle works. It works perfectly fine in the example!
• How to integrate this pressure to obtain the total force on the conductor surface (edges E3,E4,E8 & E9)?