Integrate Results of 3D FEM on mesh with arbitrary integration limits

I have results of a 3D finite element analysis stored as a list of points [x,y,z], and corresponding values V. I would like to perform 3D integration on small sub-regions of the mesh, using a piece-wise linear interpolation of this data.

I would like to be able to specify a cube in which to integrate the interpolated solution. If a portion of the cube lies outside the mesh, I would like the interpolator to return 0, instead of NaN.

I have tried the below solution, but it throws a warning (see below) and returns 0 when the integration region hangs outside the mesh. It does work on internal regions, but it is very slow (~15 seconds per integration, n_nodes ~40,000), and n_pts is ~1000. This will eventually be wrapped in lsqnonlin() for parameter estimation, so it needs to run relatively quickly for multiple iterations.

function v = RemoveNans(Si,X,Y,Z)

      v = Si(X,Y,Z);
      v(isnan(v)) = 0;

end

nodes = csvread('nodes.csv'); %3D location of points, irregularly spaced, size = [n_nodes,3]
vals = csvread('values.csv'); %solution to FEM at above nodes, size = [n_nodes,1]
int_pts =  csvread('int_points.csv'); %list of points to integrate over, size = [n_pts,3], n_pts << n_nodes

R = 1; %size of cube for integration

Si = scatteredInterpolant(nodes, vals, 'linear', 'none')

f = @(X,Y,Z)RemoveNans(Si,X,Y,Z);

integrated_values = zeros(size(int_pts));

for i=1:length(integrated_values)

      integrated_values(i) = integral3(f,...
                             int_pts(i,1)-R/2,int_pts(i,1)+R/2,...
                             int_pts(i,2)-R/2,int_pts(i,2)+R/2,...
                             int_pts(i,3)-R/2,int_pts(i,3)+R/2);

end

"Warning: Reached the maximum number of function evaluations (10000). The result passes the global error test."

Any help on getting this to work properly, or a "better" way to do this would be very helpful.