Post-process NDSolveValue result

Question

For example, after solving the Poisson equation $-\nabla^2 T = 1$ on the unit circle using NDSolveValue and the following code

eqn = -Laplacian[u[x, y], {x, y}] == 1;
sol = NDSolveValue[{eqn, DirichletCondition[u[x, y] == 0, True]}, 
  u, {x, y} \[Element] Disk[]]

how to

• sample the solution at an arbitrary point in the domain.
• integrate the "numerically obtained" temperature value over the region? This is particular relevant when solving a transient problem and one wants to know the amount of heat in the system.
• integrate the heat flux over the boundary? In other words, find the gradient normal to the local boundary and do a line integration (surface integration in 3D) over the boundary. Mathematically, the flux is $F= -\int_{\partial \Omega } \nabla T\cdot \boldsymbol {n} \,dl$, where $\boldsymbol {n}$ is the local normal direction pointing outward of the domain.

I searched for a long time without finding a solution. Thanks!

Answer 1

i) Evaluate an interpolating function:

sol[0, 0]
(* 0.25 *)

The first example of the NDSolve ref page has an example of a variant of this; many more examples are in the InterpolatingFunction ref page.

ii) Integration over the region:

mesh = sol["ElementMesh"];
NIntegrate[sol[x, y], {x, y} ∈ mesh]

(* 0.392702 *) 

iii) You could use this:

Needs["NDSolve`FEM`"]
NIntegrate[
 Grad[sol[x, y], {x, y}] . 
  BoundaryUnitNormal[x, y], {x, y} ∈ ToBoundaryMesh[mesh]]
(* -3.13394 *)

If you want to integrate only over a part of the boundary, just create a boundary mesh for that part. See the ToBoundaryMesh ref page or the ElementMesh generation tutorial. One last note, you can use things like

NIntegrate[
 NeumannValue[1, x >= 1/2], {x, y} ∈ ToBoundaryMesh[mesh]]

Let me give an example. We use this from the documentation as a starting point. There is enough text that explains the model and I am just copying it here. The interesting part will be below.

vars = {T[x, y], {x, y}};
Ω = Rectangle[{0, 0}, {0.02, 0.01}];
pars = <|"ThermalConductivity" -> 3|>;
BCTemp = 
  HeatTemperatureCondition[x == 0 || x == 0.02, vars, 
   pars, <|"SurfaceTemperature" -> 1173|>];
BCconvective = 
  HeatTransferValue[y == 0.01, vars, 
   pars, <|"AmbientTemperature" -> 323, 
    "HeatTransferCoefficient" -> 50|>];
BCradiation = 
  HeatRadiationValue[y == 0.01, vars, 
   pars, <|"AmbientTemperature" -> 323.|>];
eqn = {HeatTransferPDEComponent[vars, pars] == 
    BCconvective + BCradiation, BCTemp};
Tfun = NDSolveValue[eqn, T, {x, y} ∈ Ω];

For a plot see the documentation.

Now, we can use the boundary condition specification to compute the integrals over those parts of the boundary:

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh[Tfun["ElementMesh"]]

NIntegrate[BCconvective /. T -> Tfun, {x, y} ∈ bmesh]
(* -681.075 *)

NIntegrate[BCradiation /. T -> Tfun, {x, y} ∈ bmesh]
(*-1162.59 *)

I am not 100% sure, but I think this is what you are looking for. Your flux F is exactly the NeumannValue.