Express the double integral in terms of the function values at discrete nodes

Question

Consider the following integral:

$$ \mathcal{I}[p_{1},t] = \int \limits_{0}^{p_{1}}dp_{2}\int \limits_{0}^{p_{2}}dp_{3} \ \text{integrand} $$

where

integrand[GF_,p1_,p2_,p3_,t_]=(2 GF^2 p3^3 (-15 p3 (p1+p2)+10 (p1+p2)^2+6 p3^2) ((1-f[p1,t]) (1-f[p2,t]) f[p3,t] f[p1+p2-p3,t]-f[p1,t] f[p2,t] (1-f[p3,t]) (1-f[p1+p2-p3,t])))/(15 \[Pi]^3 p1^2);

Here, f is some unknown function which obeys the following integrodifferential equation: $$ \frac{df}{dt} - p g(t)\frac{df}{dp} = \mathcal{I}[p,t], $$ where $g$ satisfies $$ g = G \sqrt{\int p^{3}f dp} $$ Finally, $F,G$ are some constants. In principle, the second summand in the lhs of this equation may be eliminated by the change of variable $p \to k = pa$, where $a$ is defined as $$ g \equiv \frac{1}{a}\frac{\partial a}{\partial t} $$ which would also allow the solution to stay within the same domain of k at any moment of time. The system of equations represents a small piece of the equations from this question.

As far as I understand, there are no built-in methods that may solve such integrodifferential equations without the preliminary discretization (or is the method of lines applicable here?). Therefore, I need to present the equation and, in particular, the integral $\mathcal{I}$ in terms of the values of $f$ at some nodes in $p$ and then solve for the grid; they are defined as f[i,t]. The problem is that $\mathcal{I}$ is two-dimensional while the nodes are dense:

pnodes = Table[0+i*0.5,{i,1,500,1}];

Is there a command in Mathematica that automatically does this discrete representation of the integral? If not, then how to define it in order to quickly call when solving the equation?

The straightforward way may be to calculate Length[pnodes] expressions for the integral, depending on the value of p1, then export them, and then import and use them for some particular parameters and modifications entering the differential equation. My implementation is

{indexval[1], indexval[2], indexval[3]} = {i, j, k};
integranddistretized[i_, j_, k_] = 
 integrand[GF, p1, p2, p3, t]/. 
    Table[f[Symbol["p" <> ToString[l]], t] -> f[indexval[l], t], {l, 
      1, 3, 1}] /. {f[p1 + p2 - p3, t] -> f[i + j - k, t]} /. 
  Table[Symbol["p" <> ToString[f]] -> pdiscr[indexval[f]], {f, 1, 3, 
    1}]
Do[integraldiscr[i] = 
   Sum[integranddistretized[i, j, k], {j, 0, i, 1}, {k,
      0, j, 1}]
  , {i, 1, 500, 1}];

However, it is ugly and very slow, which is obvious since there are many summands for large i; also, to store this symbolic representation of the integral, I would need a lot of space: integrand[500] takes ~12 megabytes. It seems that the given realization of the approach is meaningless then.

I am wondering if there is a more elegant and faster solution.