# Controlling quality of discretized region meshes

EDIT:

Impact of this problem has been considerably lessened in v10.2. Although new version is still not performing analytical discretization at sharp corners (there are jaggies), evenness of sampling seems better and the MaxCellMeasure option is actually obeyed in the expected manner. This makes discretization much more usable.

For an example:

BoundaryDiscretizeRegion[RegionDifference[Ball[], Ball[{0, 0, 1}]],
MaxCellMeasure -> {"Area" -> 0.0025}]


I want to discretize 3D constructive solids, supported by Mathematica 10 regions functionality. Quality of resulting meshes must be controllable, and meshes must be useful for further computation.

This should be easy, but the trivial example below stuns me.

(* Simple derived region in 3D with a sharp edge: *)
r = RegionDifference[Ball[{0, 0, 0}], Ball[{0, 0, 1}]];

(* RegionPlot3D with PlotPoints handles this reasonably well: *)
RegionPlot3D[r, PlotPoints -> 50]


(* This should create a discretized mesh, and it does, but with poor quality: *)
DiscretizeRegion[r]


Edges are noticeably jaggy. This is understandable, and usually option like PlotPoints fixes these kind of issues. DiscretizeRegion has numerous options which should affect quality of its' output. I think I've played with all of them.

Biggest difference I've seen has been is through use of MaxCellMeasure, but interestingly enough, same ugly jaggies stay there even if complexity of the mesh itself increases. I would need an option like PlotPoints which apparently affects grid resolution of a marching cubes style algorithm in RegionPlot3D, but there isn't one for DiscretizeRegion. Any suggestions?

• For those looking at this question for visualization purposes I suggest taking a look at mathematica.stackexchange.com/a/48530/3056 - which, combining with RegionMember can produce very aesthetic surface renditions with nicely detailed sharp edges. – kirma Nov 2 '15 at 16:08

In principal you should be able to do

r = RegionDifference[Ball[{0, 0, 0}], Ball[{0, 0, 1}]];
rp = RegionPlot3D[r, PlotPoints -> 50];
DiscretizeGraphics[rp]


Unfortunately, this does not work and is hopefully improved in a future version.

One thing you can do, however, is use the finite element mesh generator for this:

Needs["NDSolveFEM"]
m = ToElementMesh[r,
"BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> 50},
"MeshOrder" -> 1]


You can then convert the ElementMesh to a MeshRegion simply by using

MeshRegion[m]


Note that in this case MeshRegion does not show it's output since it's beyond a threshold (too big). It has some 300000 tetrahedron elements.

You can force it to display with

Show[MeshRegion[m]]


The "MeshOrder" is set to one (= linear elements) and the "SamplePoints" correspond to the PlotPoints in RegionPlot. There is more information about boundary mesh generation in ToBoundaryMesh, ToElementMesh, ElementMesh and the related mesh generation and element mesh visualization tutorial.

One last thing I should add, is that the 3D stuff does have some rough edges.

Update version 12.1:

In version 12.1 you can use the OpenCascadeLink to do the following:

Needs["OpenCascadeLink"]
Needs["NDSolveFEM"]
r = RegionDifference[Ball[{0, 0, 0}], Ball[{0, 0, 1}]];

bmeshOC["Wireframe"[
"MeshElementStyle" -> Directive[FaceForm[Green], EdgeForm[Red]]]]


• Whoa, thanks! I'm a bit surprised FEM exposes this, instead of DiscretizeRegion, where I'd expect it to live now. Hoping to see those rough edges (pun intended) smoothed out in the future! – kirma Jul 15 '14 at 16:00
• @kirma, DiscretizeRegion should give proper results via AccuraryGoal and PrecidionGoal but that's not implemented fully for 3D. The FEM stuff exposes this to allow maximum 'hackability' and in this case that allows for a workaround. If you do find stuff that does not work, please, by all means send it to support such that it can be looked at and get fixed. Thanks. – user21 Jul 15 '14 at 16:09
• good to see you back alive, user21 :-). Is it possible to refine the grid only near the edge, in one way or another? – Albert Retey Jul 15 '14 at 18:47
• @AlbertRetey, well yes back to get the geometry and FEM stuff started ;-). You could use "MaxBoundaryCellMeasure" to refine the boundary; this works well in 2D. However, in 3D this does not yet influence the sampling of the boundary. Also, since this region has some thin parts to it it would generate quite a lot more elements. – user21 Jul 16 '14 at 6:59

Here is a way to do what you want using DiscretizeGraphics as suggested by @user21. No need to call NDSolveFEM

r = RegionDifference[Ball[{0, 0, 0}], Ball[{0, 0, 1}]]

rp = RegionPlot3D[r, PlotPoints -> 50]


Now we discretize the Graphics object with the following clever replacements that somehow works:

DiscretizeGraphics[Normal[rp /. {(PlotRange -> _) :>
PlotRange -> All, (Lighting -> _) :> Lighting -> Automatic}]]


To get even finer mesh, increase PlotPoints in RegionPlot (this will increase the time it takes to render the RegionPlot), then proceed as before. I used PlotPoints -> 100 for the image below:

• @user21, I can't seem to get the mesh as fine as that produced from NDSolveFEM in your answer. I've tried various options: MeshQualityGoal, MaxCellMeasure etc. They don't seem to have any effect. Unless I increase PlotPoints in RegionPlot which takes a while to render. Any suggestions? – RunnyKine Jul 22 '14 at 21:16
• but using NDSolveFEM is the best part of it.... ;-) I think that MaxCellMeasure is not working is a bug. MeshQualityGoal` has no much to do with the number of tets generated. It tries to generate tets that have certain quality see here – user21 Jul 23 '14 at 7:11