Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]

enter image description here

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

enter image description here

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?

share|improve this question

2 Answers 2

up vote 12 down vote accepted

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["NDSolve`FEM`"]
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]]

enter image description here

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.

share|improve this answer
    
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 at 16:00
3  
@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 at 16:09
2  
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 at 18:47
2  
@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 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]

Mathematica graphics

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

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

Mathematica graphics

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:

Mathematica graphics

share|improve this answer
    
@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 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 at 7:11

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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