NIntegrate over regions, FiniteElement options

Question

I am using NIntegrate to integrate a 3-dimensional region. When I don't set any integration options, I get a warning

The global error of the strategy GlobalAdaptive has increased more than 2000 times ... Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration

In a similar question, Can you use "MaxErrorIncreases" when numerically integrating over a region? Michael E2 suggests using the "FiniteElement" method. This integrates nicely, or at least there are no warnings or errors.

The Finite Element User Guide is very detailed but it seems to say nothing about using finite elements with NIntegrate. Neither the NIntegrate help page or the Advanced Numerical Integration guide seem to mention the method FiniteElement. Maybe it's there but I didn't find it. I would like to explore refining the mesh to check convergence but I can't find any information about how to give NIntegrate options for the FiniteElement method.

Can anyone provide or point me to some guidance on using a finite element mesh with NIntegrate?

Edit:Looking some more at the Finite Element User Guide, I see that there is a section on how to pass options for ElementMesh creation to NDSolve. I tried the same syntax for NIntegrate, setting MaxCellMeasure, and Mathematica doesn't complain, and the answer changes a bit and then converges, so maybe that works. Any confirmation would be appreciated.

Edit2: Thanks for your comments. The function I want to integrate is

f[x_, y_, z_] = (x^4 y - 2 x^3 x2 y + y^5 - 2 y^4 y2 + y^3 y2^2 + 
11 y^3 z^2 - 14 y^2 y2 z^2 + y y2^2 z^2 + 10 y z^4 - 3 y2 z^4 + 
x2^2 (y^3 - 8 y z^2) - 8 y^3 z z2 + 12 y^2 y2 z z2 - 
17 y z^3 z2 + 3 y2 z^3 z2 - 2 y^3 z2^2 + 7 y z^2 z2^2 + 
x x2 (-2 y^3 + 9 y2 z^2 + y z (-11 z + 18 z2)) + 
x^2 (x2^2 y + 2 y^3 - 2 y^2 y2 - 3 y2 z (z + 2 z2) + 
   y (y2^2 + 11 z^2 - 8 z z2 - 2 z2^2))) /(32 \[Pi]^3 (x^2 + y^2 +
   z^2)^4 ((x - x2)^2 + (y - y2)^2 + (z - z2)^2)^(5/2))

where x2, y2, z2 are O(1) parameters and I want to integrate over many values of the parameters. The region is

region = RegionDifference[RegionDifference[FullRegion[3], Ball[{0, 0, 0},ε]], 
Ball[{x2, y2, z2}, ε]]

I am interested in the limit ε << Norm[x2vector], where x2vector = {x2, y2, z2}. The integration command, without any options is

NIntegrate[f[x, y, z], {x, y, z} ∊ region]

with x2, y2, z2, and ε replaced with numbers.

Answer 1

I do not know how to solve your integral, since you did not share it. But I can fill the documentation gap about how the finite element based integration in NIntegrate works.

Here is a very rough sketch of what NIntegrate does. NIntegrate works by setting up the "LoadCoefficients" coefficient of InitializePDECoefficients to integrate over a region:

Needs["NDSolve`FEM`"]
sregion = Disk[];
exact = Pi;
vars = {x, y};
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, 
     vars}];
sd = NDSolve`SolutionData["Space" -> ToNumericalRegion[sregion]];
cdata = InitializePDECoefficients[vd, sd, "LoadCoefficients" -> {{1}}];
mdata = InitializePDEMethodData[vd, sd];
ddata = DiscretizePDE[cdata, mdata, sd];
i1 = Total[ddata["LoadVector"], 2, Method -> "CompensatedSummation"];
exact - i1
(* 2.00118*10^-6 *)

Now, if you compare that to:

i2 = NIntegrate[1, vars \[Element] sregion];
exact - i2
(* 2.60736*10^-9 *)

While the first result is already quite acceptable how can NIntegrate get a better result? The improvement is the result of adaptive mesh refinement. Adaptive mesh refinement is used for symbolic regions. That means if an ElementMesh is given to NIntegrate no adaptive mesh refinement is used; and that has manly to do with the refinement of the boundary. You can check this by comparing the low level code above with giving NIntegrate a mesh.

mesh = NDSolve`SolutionDataComponent[sd, "Space"]["ElementMesh"];
i3 = NIntegrate[1, Element[vars, mesh]];
i3 - i1
(* 0. *)

It's also roughly the same if a symbolic region is given but adaptive refinement is switched off:

i4 = NIntegrate[1, Element[vars, sregion], 
   Method -> {"FiniteElement"}, MaxRecursion -> 0];
i1 - i4

 (* -4.44089*10^-16 *)

In the adaptive case the actual integration then works very roughly in the following manner: The integrand is integrated (as with the low level FEM code above) once with "IntegrationOrder"->2 and once with "IntegrationOrder"->4. From this a per element and a total difference are computed. The integration is assumed to have converged when the iteration count is larger than 1 and and Abs of the difference between the previous result and the current result is less then a tolerance. The tolerance is computed by atol + Abs[result]*rtol, where atol is derived from AccuracyGoal and rtol from PrecisionGoal.

In case of no convergence the per element difference is used as a basis to select the elements for refinement.

To specify (mesh) options you can use:

Method -> {"FiniteElement", "MeshOptions" -> {....}}

I am not sure how useful this is going to be, since the adaptive mesh refinement is somewhat complicated and as special cases to deal with, like that boundary improvements are not automatically done in 3D, exlusion handling etc, etc..