Approximating Integral under 2D Interpolation function

Question

First, I have already been working on this problem, and have asked a similar question previously here.

I am trying to find the integrals of 2D interpolation functions that are output from NDEigensystem[]. I will be doing this integral hundreds if not thousands of times, so it is necessary to speed it up as much as possible. I do not care if the value is a scalar multiple away from the most exact value produced by NIntegrate[] and it can be (slightly) off from the a scalar multiple of the most exact value. Essentially, over many different functions the more exact value and the quick approximation should have a roughly linear relationship.

Right now the fastest NIntegrate[] function that I'm using is

fastint[eigenfunction_]:=Abs[NIntegrate[
  eigenfunction[x,y],{x,y}\[Element]eigenfunction["ElementMesh"],
  Method->{Automatic,"SymbolicProcessing"->0}]];

And the faster approximation I am using is just a sum over the "ValuesOnGrid" for the interpolation function

fastapprox[eigenfunction_]:=Abs[Total[eigenfunction["ValuesOnGrid"]]];

These two functions form a roughly linear relationship when the eigenfunctions are determined by a highly symmetric region (e.g. Disk[]). But when I do this integral with an asymmetric region, I obtain a plot similar to:

which is obviously not linear.

I feel like I am forgetting to incorporate something, but I do not know what. Again, the end goal is to either speed up NIntegrate[] or develop an approximation method that is faster while maintaining the linear relationship with a more exact value.

Any help is much appreciated, thank you.

Edit: I have developed code for a basic Reimannian approximation that is much faster than NIntegrate[] and still is able to get the results I am hoping for:

reimanapprox[region_,numberdiv_:50][functionlist_]:=
 Module[{Dx, Dy, area, squaregrid, i, inside, heights, xmin, xmax, ymax, ymin, coord},
  coord = region["Coordinates"];
  xmax = Max[coord[[All, 1]]];
  xmin = Min[coord[[All, 1]]];
  ymax = Max[coord[[All, 2]]];
  ymin = Min[coord[[All, 2]]];
  Dx = (xmax-xmin)/numberdiv;
  Dy = (ymax-ymin)/numberdiv;
  area = Dx*Dy;
  squaregrid = Flatten[Table[{x, y},
   {x, xmin, xmax, Dx},
   {y, ymin, ymax, Dy}], 1];
  inside = {};
  For[i = 1, i <= Length[squaregrid], i++,
   If[RegionMember[region, squaregrid[[i]]],
    AppendTo[inside, squaregrid[[i]]],
    False]];
  heights = Table[
   Map[eigenfunctions[[i]][#[[1]], #[[2]]] &, inside],
   {i, Length[eigenfunctions]}];
  area*Map[Total, heights]]

Which then can do many of these integrals much faster than NIntegrate[] (efunc is just a list of (30) interpolation functions that I am testing which are all over the region eregion).

RepeatedTiming[fastapprox[eregion][efunc]

(*Returns: {1.00, {-0.86141, -0.129498, -0.194252, 0.357529, 0.247833, 0.0704943, 0.0505611, ..., 0.00346365}}*)

RepeatedTiming[Table[fastintegral[efunc[[i]]], {i, Length[efunc]}]]

(*Returns: {21.1, {-0.861569, -0.129423, -0.194507, 0.357675, 0.248136,0.0706918, 0.050766,..., 0.00274389}}*)

Answer 1

Well, since an example NDEigensystem was not provided, I went ahead an took one from the docs.

{vals, funs} = NDEigensystem[
   {-Laplacian[u[x, y], {x, y}], 
    DirichletCondition[u[x, y] == 0, True]},
   u,
   {x, y} ∈ RegionDifference[Cuboid[{-3, -3/2}, {0, 3/2}], Disk[]],
   10,
   Method -> {"SpatialDiscretization" ->
     {"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> {"Area" -> 0.005}}}}}];

We can compare an "accurate" NIntegrate[] with a quadratic-order triangle rule. Warning: The triangle rule assumes the quadratic points of the element mesh are the midpoints of the sides of the triangles. This may not be true everywhere -- on a curved boundary, especially. In the present example, it can be seen not to introduce a large error.

(accurate = Table[
    NIntegrate[
     eigenfunction[x, y], {x, y} ∈ eigenfunction["ElementMesh"], 
     Method -> {Automatic, "SymbolicProcessing" -> 0}],
    {eigenfunction, funs}]; "accurate") // RepeatedTiming

Needs["NDSolve`FEM`"];
emesh = funs[[1]]["ElementMesh"];
indices = emesh["MeshElements"] /.
  {TriangleElement[t_]} :> Flatten@t[[All, 4 ;; 6]]; (* midpoints of sides *)
(fast = Table[
    Total[
     Total[Partition[eigenfunction["ValuesOnGrid"][[indices]], 3], {2}] *
      First@eigenfunction["ElementMesh"]["MeshElementMeasure"] / 3
     ],
    {eigenfunction, funs}]; "fast") // AbsoluteTiming
(*
  {2.23, "accurate"}
  {0.011379, "fast"}
*)

Okay, so it's fast. Let's now examine its accuracy.

fast
accurate
(*
  {2.08897, -3.87294*10^-7, -0.310354, 0.749835, -0.0000232558,
   0.47398, 6.61989*10^-7, -4.84232*10^-6, -0.679196, -0.290395}

  {2.08898, -3.91154*10^-7, -0.310374, 0.749833, -0.0000236342, 
   0.473997, 2.35524*10^-6, -3.86783*10^-6, -0.679223, -0.290408}
*)

NDSolve`ScaledVectorNorm[Infinity, {10^-4, 10^-5}][fast - accurate, accurate]
(*  0.483535  *)

NDSolve`ScaledVectorNorm returns an error scaled according the the {precision_goal, accuracy_goal} to be less than 1 if the error fast - accurate meets the goals, as defined by Mathematica. In short, the integrals above meet the goals corresponding to PrecisionGoal -> 4 and AccuracyGoal -> 5.