Mesh Refinement for interpolation of f[x,y] fails for all but the simplest f

Question

I would like to interpolate a function, f[x,y,z], on a mesh of the unit cube, starting with a coarse tetrahedral mesh and refining tetrahedra as needed to reduce the interpolation error below a tolerance. The problem: DiscretizeRegion[Cuboid[],MeshRefinementFunction->mrf] ignores my mrf. In an attempt to understand, I have simplified to a toy version of the problem (2D, on a square instead of a cube, with a simple form for f). I discovered that of the following 5 mathematically equivalent forms for f, the last works and the others do not:

fDirect[{x_, y_}] = Cos[2*Pi*x]*Sin[2*Pi*y]
fDelayed[{x_, y_}] := Cos[2*Pi*x]*Sin[2*Pi*y]
fEmbedded1 = Function[{p}, fDirect[p]]
fEmbedded2 = fDirect[#1] & 
fPureFun = Function[{p}, Cos[2*Pi*p[[1]]]*Sin[2*Pi*p[[2]]]]

Here is what the unrefined mesh looks like:

sq = DiscretizeRegion[Rectangle[], MaxCellMeasure -> Infinity]

The resulting mesh has 4 nodes and 2 triangular elements.

Now we define a refinement function. For the toy problem I have adopted a simplified form of the error estimate. I check the error only at the center of each tetrahedron. The center is at Mean[vertices]. The true value there is f[Mean[vertices]]. My "interpolated value" there is just the mean of f at the individual vertices, Mean[f/@vertices]. When Abs[error]>tolerance, the function returns True. A diagnostic print statement lets us see when the mesh refinement function is called.

This does not work for the first definition of f. Here is what that looks like:

 Print["fDirect"]
 mrf = Function[{vertices, area}, 
   Module[{fest, ftrue, error, tolerance, test}, 
    fest = Mean[fDirect /@ vertices]; 
    ftrue = fDirect[Mean[vertices]]; 
    error = fest - ftrue; tolerance = 0.43; 
    test = Abs[error] > tolerance; 
    Print[{vertices, area, fest, ftrue, error, tolerance, test}]; 
    test
    ]
   ]; 
sq = DiscretizeRegion[Rectangle[], MeshRefinementFunction -> mrf, 
  MaxCellMeasure -> Infinity]
Print["The mesh has ", Length[mc = MeshCoordinates[sq]], " nodes."]
Print["The mesh has ", 
 Length[mcells = MeshCells[sq, 2]], " triangular cells."]

Notice that the mrf is called only once (there is output for only one diagnostic Print[]) despite that there are two triangles. In that one call the error exceeds tolerance and the mrf returns true (i.e., do a refinement), yet there is no subsequent refinement. Definitions 2 - 4 produce the same result.

It does work for definition 5. Here is what it's supposed to look like when it works:

mrf = Function[{vertices, area}, 
   Module[{fest, ftrue, error, tolerance, test},
    fest = Mean[fPureFun /@ vertices];
    ftrue = fPureFun[Mean[vertices]];
    error = fest - ftrue;
    tolerance = 0.43;
    test = Abs[error] > tolerance;
    Print[{vertices, area, fest, ftrue, error, tolerance, test}];
    test]];
sq = DiscretizeRegion[Rectangle[], MeshRefinementFunction -> mrf, 
  MaxCellMeasure -> \[Infinity]]
Print["The mesh has ", Length[mc = MeshCoordinates[sq]], " nodes."]
Print["The mesh has ", 
 Length[mcells = MeshCells[sq, 2]], " triangular cells."]

In the working case, we get a diagnostic print output for each cell of the mesh with repeats for several iterations, until all triangular elements have interpolation error estimates below tolerance.

It seems odd that the change that makes a difference here is a difference in the definition of f, which is internal to a module within a function definition. This seems a failure of encapsulation.

Definitions 3 and 4 were attempts to hide one of the apparently unacceptable forms of definition inside of the acceptable pure function definition form, but this did not work.

Unfortunately, my non-toy version of f is not easily written as a pure function (the one mode that worked above). That function is the end result of tens of pages of development within the notebook. It calls a number of functions, uses a numerical differential equation solver, an iterative root finder,... It's not a simple function, which is why I want to tabulate and interpolate it.

Is this a bug? Is there a workaround--another way to access mesh refinement within Mathematica? Or would this work if only I did it correctly? (How does one do it correctly?)

Answer 1

You could precompile the mesh refine function with additional information.

For example,

cmrf = Compile[{{vertices, _Real, 2}, {area, _Real}}, 
  Module[{fest, ftrue, error, tolerance, test}, 
   fest = Mean[fDirect /@ vertices]; ftrue = fDirect[Mean[vertices]]; 
   error = fest - ftrue; tolerance = 0.43`; 
   test = Abs[error] > tolerance; test], {{fDirect[_], _Real}}]