How to speed up NIntegrate when using an interpolating function output from NDEigensystem?

Question

I have been working with NDEigensystem in order to find the resonant frequencies of different shaped drums. I would also like the coefficients of the different eigenfunctions so that each of these frequencies has a corresponding intensity (the coefficients squared). The only problem is that the NIntegrate command with my two interpolating functions is running really slow. First, here is my function to model a struck drum:

pluck[{{f_, g_}, {umin_, umax_}}, Nbordpts_][{Ax_, Ay_}] := 
  Module[{hmax, d, l, recurt, constructxy, oneline, datapts, recurtpoints, mesh},
  hmax = 1;
  d = .05;
  (* l:[0,1]->R^2 *)
  l[{x1_, y1_}, {x2_, y2_}][t_] := {(1 - t) x1 + t x2, (1 - t) y1 + t y2};
  (* Constructs points for a line from (0,0) to (1,hmax) *)
  recurt[{n_, t_, z_}] := {n + 1,
   1 - (1 - t)/(hmax + d/(n + 1) - z) (2 d)/(n + 1),
   hmax - d/(n + 1)};
  recurtpoints = Join[NestList[recurt, {0, 0, 0}, 10],
   Table[
   {0, (1 - (2 d)/(hmax + d)) i/4, (hmax + d) (1 - (2 d)/(hmax + d)) i/4},
   {i, 1, 3}]];
  (* constructxy: constructs points in R^3 from points recursively determined previously
  in recurt using a border point calculated from coordinate function heads f,g *)
  constructxy[U_][{n_, t_, z_}] := {l[{f[U], g[U]}, {ax, ay}][t], -z};
  (* oneline: function to create set of points {{x,y},z} along one line from
  {f[U],g[U]} to {ax,ay} - f,g need to be the *heads* of the coordinate functions
  e.g Cos,Sin NOT Cos[u],Sin[u] *)
  oneline[Num_][U_] := (Map[constructxy[U], recurtpoints]);
  (* Essentially create these interpolation points for a bunch of lines from a
  bunch of different border points and put them all into one set, then
  interpolate hmax is Absolute value of maximum depth (default is 1)*)
  datapts = (DeleteDuplicatesBy[Append[Flatten[Map[oneline[10],
   Range[umin, umax, (umax - umin)/Nbordpts]], 1], {{ax, ay}, -hmax}],
   First] /. {ax -> Ax, ay -> Ay});
  (* This requires a random addition of a small amount to make sure that
  Mathematica does not mistake a triangle in the mesh as having zero area *)
  mesh = NDSolve`FEM`ToElementMesh[datapts[[All, 1]] + RandomReal[{-10^-9, 10^-9},
   {Length[datapts]}]];
  NDSolve`FEM`ElementMeshInterpolation[{mesh}, datapts[[All, 2]], InterpolationOrder -> 1]];

This function creates a list of points in the form {{x,y},z} and then constructs a first order interpolation of them.

{circlevals,circlefuncs}=NDEigensystem[{-Laplacian[f[x, y], {x, y}],
 DirichletCondition[f[x, y] == 0, True]}, {f}, {x, y} \[Element]Disk[], 40];
circlefreqs=Sqrt[circlevals];
circlepluck=pluck[{{Cos,Sin},{0,2 Pi}},100][{0,0}];

Now the list {circlefreqs,circlefuncs} contains the eigenfrequencies and the eigenfunctions for the region of the unit disk, and circlepluck contains an interpolating function that models a circular drum struck at its center. Here is the command that consistently takes over two minutes (everything else is fast on its own):

Table[
 NIntegrate[circlepluck[x, y]*circlefuncs[[i]][x, y], {x, y} \[Element]Disk[],
 {i, 1, Length[circlefuncs]}]

I have tried changing the interpolation order in NDEigensystem as well as many of its other options such as making the discretizing of the region more precise and turning off VectorNormalization. Also, I have tried different methods for NIntegrate (like "Trapezoidal", "MonteCarlo", among others) and changing the precision/accuracy goals, but all of them are either too far away from the actual values or too slow.

I do not need the values to be incredibly accurate, but I would like their relative ratios to stay somewhat accurate (To explain slightly more broadly: I will be running this code for many different borders and I would like to find out which border sounds the best).

Any help would be much appreciated, thank you.

Answer 1

Making my comment an answer:

Try to extract the mesh from the eigenfunction and use that for integration. Something like:

NIntegrate[ circlepluck[x, y]*circlefuncs[[i]][x, y], {x, y} \[Element] circlefuncs[[1]]["ElementMesh"]]

This would switch off the adaptive mesh refinement.