Integration of solution in different subdomains

Hello,

I know that this topic has been somehow discussed in other threads, but it didn't fully answer my questions. I want to find the magnetic energy stored in subdomains in a magnetostatic solution. That is why I need to integrate over the absolute of the gradients of the solution in a specific subdomain. I didnt want to use the meshgrid function and evaluategradient since it would introduce another error in the solution. That's why I thought I could use the value at of the magnetic flux at the midpoints of the triangles and then multiply those with the area of the triangle for all triangles in a subdomain. This worked somehow but proposed another error since the triangles that are next to another subdomain have one or two nodes that also belong to another subdomain. That results in a wrong solution, especially for small subdomains.

My question:

1. Since the solver finds the solution for the nodes, how does it choose if the node is in one or the other subdomain when the node is on the border between those subdomains.

2. How can I exclude those nodes or is there a better way to integrate over a subdomain.

My Code:

[p,~,t] = meshToPet(mesh);
%Get the right area for integration for all named areas in sdl
tsub = t(:,t(4,:)==sdl(1));
%Get the interpolated midpoint
ut = pdeintrp(p,tsub,u);
[a,~,~,~] = pdetrg(p,tsub);
a_sum=sum(a);
%Get mean value of solution in the area
u_mean = sum(ut.*a)/a_sum;

In my case the solution u is the absolute of the gradient of the NodalSolution.

Thank you in advance!