# Integrate result of finite element calculation in MATLAB

I have a PDE that I solve in MATLAB over a 2D domain using the "recommended workflow", i.e. the functions descg, createpde and so on. The solver (solvepde) works and return a solution which looks reasonable. So far so good.

Now I want to take the result and calculate various integrals on it over the whole domain. Mainly, I would like to calculate $$\int_\Omega \Big|\vec\nabla u\Big|^2 d^2x$$ where $\Omega$ is the 2D PDE domain and $u$ is the solution of the PDE.

What would be the best way to achieve this?

---- More details ----

If it makes any difference, the equation is $$\nabla^2 u -a u=0$$ with Neumann boundary conditions. $a$ is some small number that I put in just to make the solution unique (since all my boundary conditions are Neumann the solution is defined only up to an additive constant. If you have a better way to solve the non-uniqueness then I'll be glad to hear about it too). The domain is a circle with two circular holes. Here's the code that generates the geometry and the solution:

gd=[1 0 0 1; 1 .8 0 .1; 1 -.4 0 .1]';
ns=char('bound','p1','p2')';
sf='(bound-p1)-p2';
[dl, bt]=decsg(gd,sf,ns);
model=createpde();
geometryFromEdges(model,dl);
applyBoundaryCondition(model,'neumann','Edge',1:12,'q',0,'g',1);
specifyCoefficients(model,'m',0,'d',0,'a',1e-6,'f',0,'c',1);
generateMesh(model,'Hmax',.06);
u=model.solvepde;


If the domain $\Omega$ is nice, the solution to the pure Neumann problem is defined up to an additive constant. And any conjugated gradient method for solving $Ax=b$ with such singular A will converge if you make the righthand side orthogonal to the kernel which is the constant vector (i.e., subtract from b the average value of its elements).
So the only thing to do then is to make the solution obtained from the linear solver also orthogonal to constant vector. For me this way seems to bit a little bit simpler than adding small perturbation to the problem via $au$.
I am sure that matlab either has a function in this toolbox, to compute this term (it is part of the $H^1$ norm) or a function that will give you the stiffness matrix A of your FEM discretization. Then you compute this integral as u'*A*u.