NIntegrate InterpolatingFunction Kernel Crash

Question

I would like to Investigate if the edges of this 2D heat equation are insulated by the Mathematica assumed Neumann conditions.

aap = NDSolve[{\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(f[x, y, t]\)\) == 0.05 \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(f[x, y, t]\)\), 
    f[x, y, 0] == E^(-40 (x^2 + y^2))}, 
   f, {t, 0, 1}, {x, y} \[Element] Disk[]][[1]]

The problem is that Mathematica seems unable to Integrate the solution to the problem.

    Export["test.gif", 
 Table[Plot3D[aap[x, y, i], {x, y} \[Element] Disk[], 
   PlotRange -> {Automatic, Automatic, {0, 2}}], {i, .1, 5, 0.2}], 
 "AnimationRepetitions" -> \[Infinity], "ImageSize" -> 20]

I try

NIntegrate[f[x, y, 0] /. aap, {x, y} \[Element] Disk[]]

But my kernel just quits, same thing on Mathematica Online.

Do you also see this?/Is someone able to help?

edit:

I now did get the result one time:

In[5]:= NIntegrate[f[x, y, 0] /. aap, {x, y} \[Element] Disk[]]

During evaluation of In[5]:= Outer::heads: Heads Alternatives and List at positions 3 and 2 are expected to be the same.

During evaluation of In[5]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

During evaluation of In[5]:= NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.07854482887754174` and 2.6180493354023447`*^-7 for the integral and error estimates.

Out[5]= 0.0785448

which is the result I need. After that I managed to crash the kernel again and I am back at the same point...

edit (2) the workaround solution as mention below works! The following now shows there is indeed no outflow of "heat".

aap = NDSolveValue[{D[f[x, y, t], t] == 
     0.05 Laplacian[f[x, y, t], {x, y}], 
    f[x, y, 0] == 40/\[Pi] E^(-40*(x^2 + y^2))}, 
   f, {t, 0, 10}, {x, y} \[Element] Disk[]];

and plot

Plot[NIntegrate[
  aap[x, y, t], {x, y} \[Element] aap["ElementMesh"]], {t, 0, 10}, 
 PlotRange -> {0, 2}]

to get:

Answer 1

That's a bug which I filed. As a workaround you can use:

aap = NDSolveValue[{D[f[x, y, t], t] == 
     0.05 Laplacian[f[x, y, t], {x, y}], 
    f[x, y, 0] == E^(-40*(x^2 + y^2))}, 
   f, {t, 0, 1}, {x, y} \[Element] Disk[]];
NIntegrate[aap[x, y, 0], {x, y} \[Element] aap["ElementMesh"]]
0.07854480753557473`