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;