I have a PDE that I solve in
MATLAB over a 2D domain using the "recommended workflow", i.e. the functions
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;