Expensive MeshRefinementFunction

I want to discretize a region using DiscretizeRegion. Now, I want the mesh to be more refined, when the curvature of a certain function f is larger than some fixed value, e.g.

DiscretizeRegion[
Interval[{0, 1}],
MeshRefinementFunction -> Function[{v, len},
(* sample more points when curvature is large *)
Abs@First@Differences[f /@ {v[[1]], Mean[v], v[[2]]}, 2] > 0.005
]
]


The problem is, that f is very expensive to evaluate. Therefore, I want to store all the generated values, so I can use them later. Naturally, I use the MeshRefinementFunction option, and Sow (and Reap) for this, but it doesn't work. As soon as Sow or Print appears in the function, MMA generates an error message and DiscretizeRegion does not evaluate. Here is a minimal example:

DiscretizeRegion[
Interval[{0, 1}],
MeshRefinementFunction -> (Sow@#2 > 0.1 &)
]


The messages are:

1. MeshRegion::dgcell: "The cell Line[{1,2}] is degenerate."
3. DiscretizeRegion::drf: "DiscretizeRegion was unable to discretize the region Interval[{0,1}]."

Everything is fine as soon as you remove the Sow. I suspect, that DiscretizeRegion tries to compile this function... What can I do about this (except writing my own version of DiscretizeRegion)?

• This does not seem right, I filed it as a bug. As a side note, the missing message text suggest that you are not using the latest Version (V10.3) That version has many improvements for Comp-Geom (though this issue seems to persists) – user21 Nov 7 '15 at 7:28
• Ok, thank you @user21! – JEM_Mosig Nov 7 '15 at 7:33

This does not solve the problem, but provides some additional observations.

As

DiscretizeRegion[Interval[{0, 1}], MeshRefinementFunction -> (Echo@#2 > 0.01 &)]


or

DiscretizeRegion[Interval[{0, 1}], MeshRefinementFunction -> ((Print[#2]; #2 > 0.01) &)]


work without any error messages I conclude that Reap is used somewhere within DiscretizeRegion and therefore using Sow causes these error messages.

One could use

list = {};
DiscretizeRegion[Interval[{0, 1}],
MeshRefinementFunction -> ((AppendTo[list, #2]; #2 > 0.01) &)]


instead of Sow.

However, the resulting output differs from using

DiscretizeRegion[Interval[{0, 1}], MeshRefinementFunction -> (#2 > 0.01 &)]


and the print outs show that a default 0.05 is used irrespective of the specified value.

The same behavior can be observed when using

DiscretizeRegion[Line[{{0}, {1}}],
MeshRefinementFunction -> Function[{v, a}, Print[a]; a > 0.01]]


or

DiscretizeRegion[Interval[{0, 1}], MeshRefinementFunction -> ({##}[[2]] > 0.001 &)]


An error message that reveals that Compile is being used by DiscretizeRegion can be forced with

DiscretizeRegion[Interval[{0, 1}], MeshRefinementFunction -> (##[[2]] > 0.001 &)]


These problems seem to be absent for 2D regions:

Reap@DiscretizeRegion[Disk[], MeshRefinementFunction ->
Function[{vertices, area}, Sow[area] > 0.0005 (1 + 10 Norm[Mean[vertices]])]]

• Thank you for your comment! Unfortunately, your code does not work properly on MMA 10.1: it does not produce any messages, but DiscretizeRegion now ignores the MeshRefinementFunction option completely. When you replace 0.1 by 0.01, you will see that you get different results with and without the AppendTo part. Or perhaps this is a bug in MMA 10.1 (which I am still using)? Does it work properly on MMA 10.3? – JEM_Mosig Nov 7 '15 at 7:18
• Also, it should not matter if they used Sow inside, especially when I use a unique tag, which I tried. – JEM_Mosig Nov 7 '15 at 7:20
• @user, Echo[] is a new function in 10.3, so your version won't have it. It is indeed bizarre that the thing still doesn't work even after adding a tagged Sow[]. – J. M. is computer-less Nov 7 '15 at 7:22
• I was refering to the second code. So with ((AppendTo[list, #2]; #2 > 0.02) &) DiscretizeRegion behaves as if MeshRefinementFunction was not specified, but with ((#2 > 0.02) &) it works all fine. – JEM_Mosig Nov 7 '15 at 7:27
• @J.M. Using a tag for Sow doesn't make a difference, when no tag was used for the hypothetical Reap inside DiscretizeRegion. For example, try Reap[Reap[Sow[x, "tag00"]];] and Reap[Reap[Sow[x, "tag00"], "noTag"];]. – Karsten 7. Nov 7 '15 at 7:52