0
$\begingroup$

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]

Coarse mesh of unit square without refinement

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

Failed mesh refinement using a direct definition of f[x,y]

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

Working mesh refinement with f defined as a pure function

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?)

$\endgroup$
2
  • $\begingroup$ Change these two lines to: fest = Mean[fDelayed @@@ vertices]; ftrue = fDelayed @@ Mean[vertices]; and use this modification of fDelayed instead fDelayed[x_, y_] := Cos[2*Pi*x]*Sin[2*Pi*y] and it seems to work for me. $\endgroup$ – flinty Jul 1 '20 at 15:07
  • $\begingroup$ You are right. This works for me too. Do you know why some of these forms work and others do not? It may be that my actual f[x,y,z] breaks some of the unknown rules. It would be easier to avoid this if the rules were more explicit. $\endgroup$ – John V. Jul 1 '20 at 17:57
1
$\begingroup$

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}}]
$\endgroup$
2
  • $\begingroup$ Do you say PREcompile because the mrf must be compiled and is otherwise compiled under the hood? (The documentation doesn't mention this.) Your compile idea also works for me. I will see if my actual f[x,y,z] will compile. According to the Compile[ ] documentation, if it does not compile, Mathematica will revert to ordinary Wolfram Language code. Ordinarily this should only affect speed of execution, but in this case I'm afraid it is likely to affect whether the MeshRefinementFunction option works or not. $\endgroup$ – John V. Jul 1 '20 at 18:12
  • $\begingroup$ Thanks. The explicit Compile[ ] works on my real (non-toy) f[x,y,z]. We have achieved what amounts to success: the ability to move on to the next problem. (The next problem is that it takes hours to refine my mesh even though the tolerance is set much higher than I would like.) I am upvoting and accepting this answer. $\endgroup$ – John V. Jul 3 '20 at 13:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.