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..