# Generating mesh from region - possible bug

I want to perform a calculation over a region that is an isosceles triangle with a circular segment attached to the base. The region generalises a circular sector, where the apex of the triangle is not necessarily the centre of curvature. I use:

ct = 57*Pi/180.;(* particular case *)
t0 = 20*Pi/180.;(* particular case *)
m = Tan[t0];
bigtri = Triangle[1.2*{{0., 0.}, {1., m}, {1., -m}}];
smalltri = Triangle[Cos[t0]*{{0., 0.}, {1., m}, {1., -m}}];
circ = Disk[{Cos[t0] - Sin[t0]/Tan[ct], 0.}, Sin[t0]/Sin[ct]];
region = RegionUnion[RegionIntersection[bigtri, circ], smalltri];
RegionPlot[region]


(Explanation, for completeness: t0 is half the apex angle, ct is the "contact angle" between the base and the arc, and the slope sides are unit length. The triangle is on its side, so the base is to the right. I find the intersection between the full circle and an oversize triangle to chop off those parts of the circle that are outside the required region, then I take the union of the result with the correctly-sized triangle otherwise I lose the apex in some more-curved cases.)

RegionPlot shows me the shape I was expecting:

Now to generate a mesh and take a look at it:

Needs["NDSolveFEM"];
ToElementMesh[region]["Wireframe"]


and it looks fine:

But, if I change the contact angle ct slightly from 57 to 55 degrees:

ct = 55*Pi/180.;(* particular case *)


and repeat with everything else unchanged, I get the expected RegionPlot (almost the same as before) but a defective mesh:

Try again with ct half way between, 56 degrees, and Mathematica hangs for several seconds before the kernel kills itself!

I found quite a few cases of this kind of behaviour, for unpredictable combinations of parameters.

I'm assuming this behaviour is a bug. (If not, why not?) Are there ways to reliably avoid it that don't have their own unpredictable sets of problematic parameters?

I'm running version 11.1.0.0 on 64-bit Windows 8.1 Pro.

• That's a bug in the continuation code that is fixed in the next release. – user21 Jul 25 '17 at 17:32
• OK thanks - I'll look forward to the next release. – pystab Jul 31 '17 at 14:37

ct = 55*Pi/180.;(*particular case*)t0 =  20*Pi/180.;(*particular case*)m = Tan[t0];
bigtri = Triangle[1.2*{{0., 0.}, {1., m}, {1., -m}}];
smalltri = Triangle[Cos[t0]*{{0., 0.}, {1., m}, {1., -m}}];
circ = Disk[{Cos[t0] - Sin[t0]/Tan[ct], 0.}, Sin[t0]/Sin[ct]];
region = RegionUnion[RegionIntersection[bigtri, circ], smalltri];
RegionPlot[region];

Needs["NDSolveFEM"];
ToElementMesh[region, "BoundaryMeshGenerator" -> "RegionPlot",
MaxCellMeasure -> {"Length" -> .1}]["Wireframe"]


• Could you comment on the exact changes you made in the code? It would improve readability: there seems to be a lot of overlap between your code and the OP's, so one has to hunt for change. – MarcoB Jul 10 '17 at 15:28
• @MarcoB The OP has solved with ct = 57*Pi/180 and cannot solve with ct = 55*Pi/180 (see his pictures and his remark "if I change the contact angle ct slightly from 57 to 55 degrees:"). In ToElementMesh, I added the necessary options. – user36273 Jul 10 '17 at 16:01
• The option that seems to work here is "BoundaryMeshGenerator", since removing the MaxCellMeasure doesn't change the result. Also note that your code contains a RegionPlot[region]; that doesn't seem to accomplish anything. I'd suggest trimming the unnecessary bits off, again for the sake of readability, and adding your commentary to your answer itself. – MarcoB Jul 10 '17 at 16:39

Well that's a bit odd. Someone posted an answer to my question yesterday, but it's disappeared without trace. At least I can't see it below my question today, whereas it was there yesterday. It was just some code, with no explanation or commentary.

Anyway I did make notes and, for the benefit of anyone who may find it useful, the missing answer simply replaced the last line in my example with:

ToElementMesh[region, "BoundaryMeshGenerator" -> "RegionPlot",
MaxCellMeasure -> {"Length" -> .1}]["Wireframe"]


I found the assignment to MaxCellMeasure to be unnecessary and it's "BoundaryMeshGenerator" that does the business. According to the documentation for ToBoundaryMesh, the default "Continuation" method is supposed to be better at resolving corners and cusps (not in my case!) while the "RegionPlot" method (used above) is faster but worse at corners.

I believe there is a bug with the "Continuation" method and for now I must forgo the advantages it should offer, but at least I can make progress with the workaround.

To whoever posted the missing answer: thanks, and I'm sorry I can't up-vote it. Of course if it was deleted because it was defective in some way then I'd like to know, though it looks good for now. Otherwise - why did a useful (and inoffensive) post simply vanish?

• I'm sorry, but I leave StackExchange for various reasons. I have added maxCellMeasure extra, in order to be able to refine the grid if necessary. This is my last comment. – user36273 Jul 11 '17 at 14:17