# NDSolve post-processing: Calculate the flow over a FEM-boundary

I'm trying to calculate the heat loss from buried pipes with NDSolveFEM. For this I set Dirichlet conditions at the model top and bottom and temperature dependant Neumann conditions at the contact surface between the fluid inside the pipes and the pipe. The solution looks reasonable, but I have trouble calculating the flow resulting from the Neumann conditions. The steady state solution looks like this:

Here's what I've done: This is the geometry, mesh and coefficients: http://s000.tinyupload.com/index.php?file_id=46720655068131319221

{boundaryMeshs, bmesh, mesh, κTensor, steadyState} =
Import[NotebookDirectory[] <> "problemData.mx"];


I started by assembling the boundary mesh from multiple boundary meshs. E.g. boundaryMeshs[1] being the right contact area fluid-pipe, which is the one of interest. The κTensor includes the thermal diffusivity of each component and has the unit [m^2/s].

op = Inactive[ Div][-κTensor.Inactive[Grad][T[x, y], {x, y}], {x, y}]

ΓDirichlet = {
DirichletCondition[T[x, y] == 10, ElementMarker == 9],
DirichletCondition[T[x, y] == 10, ElementMarker == 11]
};

ΓNeumann = NeumannValue[(55 - T[x, y]) 25, ElementMarker == 1] +
NeumannValue[(30 - T[x, y]) 25, ElementMarker == 5];


ElementMarker==1 corresponds to the mentioned area and the flow is defined by the fluid temperature, the pipe temperature and the heat transmisson coefficient which is calculated by the Nusselt number, but considereed to be 25 for simplification here.

steadyState = NDSolveValue[{op == ΓNeumann, ΓDirichlet}, T, {x, y} ∈ mesh]


Now I want to calculate the actual resulting flow, preferably by using my mesh. For this I picked the boundary mesh corresponding to this contact area and use the same equation as for the Neumann condition.

NIntegrate[(55 - steadyState[x, y]) 25, {x, y} ∈ boundaryMeshs[[1]]]
Out= 0.0000292133


The result would be 0.03 mW/m_pipe if I am interpreting this right. This doesn't seem right. Earlier I had the same caluclation resulting in 0 a lot of times. I guess my Integral is over the area of boundaryMeshs[1], which is a RegionHole and not over its surface. I just can't figure out how to do this right. I guess there's a very easy/elegant way. If somebody has an idea and could help me, I would be very grateful. I would expect the result to be of the magnitude 2-20 W/m.

ADDITION: If I calculate the same for the left pipe I get an even smaller number:

NIntegrate[(30 - steadyState[x, y]) 25, {x, y} \[Element]
boundaryMeshs[[5]]]
Out: 1.19173*10^-6


I calcluated the flow for the right pipe by defining my integration limits "by hand" and get a result in the range I expect, but it would be much more efficient for me to use the boundary markers, since I want to change the geometry.

NIntegrate[(55 -
`