# Laplacian or Grad of an InterpolatingFunction

For obtaining a Schlieren image from an equation for density, I need to calculate the first derivative of density and make a contour plot.

The below code snippet calculates the first derivative w.r.t x and y and makes a contour plot. All is well until this step. hSol is an interpolating function polynomial, TRup is the final time of calculation/rupture time in NDSolve and fac is just a value between 0 and 1.

T = 1 - hSol[x, y, fac*TRup]/(hSol[x, y, fac*TRup] + 1);
a = 9.2/10^4;
b = 4.5/10^7;
r = 0.97/(1 + a*T + b*T^2);
ContourPlot[Evaluate[D[r, x] + D[r, y], {x, 0, L}, {y, 0, L}],
PlotRange -> {{0, L}, {0, L}, {0, 3.5}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 18},
Contours -> 75]


If I replace the dx dy terms with Grad[r] or Laplacian[r], it errors out:

1. Either it locks up the RAM and my computer crashes.
2. Or All the code in my notebook turns to BLUE as if it were never evaluated in the first place.

I tried it with a simpler problem as in Example 1 and i could replace y' in the ParametricPlot of In[3] with Laplacian[] and things worked well. (I do realize that y' is not Laplacian...)

*(MY CODE)*Mathematica code that produces the interpolating function polynomial:

(*"Specifies the number of previous lines of input and output to keep in a Mathematica session."*)
\$HistoryLength=0;
(*Loads necessary packages to perform vector analysis*)
Needs["VectorAnalysis"]
Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
(*Clears equation symbols*)
Clear[Eq0,EvapThickFilm,h,Bo,\[Epsilon],K1,\[Delta],Bi,m,r] (*Eq0 and EvapThickFilm are variable names. These can be changed as long as the change is made through the entire document*)
(*'h' is the film thickness variable, 'Bo' is the Bond number (G/S), '\[Epsilon]' is the epsilon from the maximizing wavelength equation, 'K1' is the non-equilibrium coefficient, '\[Delta]' is the delta term from the max wavelength equation, 'Bi' is the Biot number, 'm' is the m term and 'r' is the modified Rayleigh number term*)
Eq0[h_,{Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_}]:=\!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]h\)+Div[-h^3 Bo Grad[h]+h^3 Grad[Laplacian[h]]+(\[Delta] h^3)/(Bi h+K1)^3 Grad[h]+m (h/(K1+Bi h))^2 Grad[h]]+\[Epsilon]/(Bi h+K1) + (r)D[D[(h^2/(K1+Bi h)),x] h^3,x] ==0; (*This is the evolution equation in 3D cartesian coordinates*)
SetCoordinates[Cartesian[x,y,z]];
EvapThickFilm[Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_]:=Eq0[h[x,y,t],{Bo,\[Epsilon],K1,\[Delta],Bi,m,r}];

(*Domain size in non dimensional units*)L=76.1311;  (*Maximum time of iteration*)TMax=5000*100;
Off[NDSolve::mxsst];
Clear[Kvar];
(*Kvar[t_]:=  Piecewise[{{1,t<=1},{2,t>1}}]*)
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol=h/.NDSolve[{
(*Bo,\[Epsilon],K1,\[Delta],Bi,m,r*)
(*Non dimensional number values are fed into the evolution equation which is called here with 'EvapThickFilm[...]'. The arguments, in order are, Bo, \[Epsilon], K1, \[Delta], m, r. These are the exact values used in Oron et al.*)
EvapThickFilm[0.0,0,1,0,1,0.0544911,0],

(*Periodic boundary conditions at the left and right walls in x and y directions*)
h[0,y,t]==h[L,y,t],
h[x,0,t]==h[x,L,t],
(*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)
(*Initial condition as used in Oron et al.*)
h[x,y,0]==1+(-0.05 Cos[2\[Pi] x/L] -0.05 Sin[2\[Pi] x/L])Cos[2 \[Pi] y/L]
},
h,
{x, 0, L},
{y,0, L},
{t, 0, TMax},
Method->{"BDF","MaxDifferenceOrder"->1},
MaxStepFraction->1/50
][[1]]

fac = 0.99;
hGrid = InterpolatingFunctionGrid[hSol];
{TMin, TRup} = InterpolatingFunctionDomain[hSol][[3]]
Length[hGrid]
{nX, nY, nT} = Drop[Dimensions[hGrid], -1]

-

You'll have to make sure you plot a real scalar quantity if you want to leave all the other parts of the ContourPlot unchanged. So for example you could replace D[r, x] + D[r, y] by Grad[r].{1,1,0} but not directly by Grad[r].
I see some other issues - the evaluate should wrap only the function to be plotted, not the variable range, i.e., Evaluate[D[r, x] + D[r, y], ... should be Evaluate[D[r, x] + D[r, y]]`.