Extract symbolically processed form of NIntegrated function

As I understand NIntegrate does some symbolic preprocessing of inputs before it actually goes ahead and numerically integrates. I’m wondering if there is some way to access this intermediate preconditioned form.

Let me demonstrate:

expr[a_, b_, c_] = -((a xp + b yp - c zp)/Sqrt[(0.5` - xp)^2 + yp^2 + zp^2]);
dom = Triangle[{{0.`, 0.`, 1.`}, {0.25, 0.8, 0.4}, {-0.7, 0.5, 0.4}}];
Timing[NIntegrate[expr[1,2,3], {xp, yp, zp} [Element] dom]]
Timing[NIntegrate[expr[1,2,3], {xp, yp, zp} [Element] dom, 
  Method -> {Automatic, "SymbolicProcessing" -> 0}]]

{0.03125, 0.474964}
{0.21875, 0.47448}

As you can see there’s a healthy performance gain when allowing NIntegrate to precondition its inputs.

I need to compute thousands+ of similar integrals over different domains and different parameters. It would be nice if I could get access to this preconditioned form so I could inject my parameters directly into that or perhaps compile it, rather than duplicating this simplification step each time.

NIntegrate[expr[1, 2, 3], {xp, yp, zp} [Element] dom, 
 MaxRecursion -> 0, PrecisionGoal -> 1, IntegrationMonitor -> Print]

(* During evaluation of In[35]:= {NIntegrate`IntegrationRegion[{{0,0},{1,1}},Experimental`NumericalFunction[{xp,yp},-((3.21926 (1-xp) (-0.821918+1. xp+0.575342 (1-xp) yp))/Sqrt[1.17647 -1.36471 xp+1. xp^2-0.470588 (1-xp) yp+1.10118 (1-xp) xp yp+1.03529 (1-xp)^2 yp^2]),-NumericalFunctionData-],{},NIntegrate`MultiDimensionalRule[{{{{0.,0.},{0.179284,0.},{0.474342,0.},{0.474342,0.474342},{0.344124,0.344124}},{-0.193873,0.149367,0.0518214,0.0101611,0.0871183},{1.13809,-0.354748,0.00723975,-0.0241325,0.0871183}},{{},{},{},{}},{2,5,-1}}]]} *)

(* 0.474984 *)