# What is causing these artifacts in RegionPlot3D?

I have a simple program like this:

length = 10;
o = {0, 0, 0};
dp = 0.25 (*m*);
rho = 2345 (*kg/m^3*);
rho = dp^3*rho;
region = Cuboid[o, o + length];
cavity = Table[Ball[RandomPoint[region], RandomReal[{0.6 dp, 1.5 dp}]], 4];
region = Fold[RegionDifference, region, cavity];
RegionPlot3D[region, PlotStyle -> Opacity[0.5], PlotPoints -> 80]


What is causing the artifacts that are marked with red arrows? I think the problem is derived from PlotPoints, but how can I avoid this?

## 3 Answers

It looks like a bug in RegionPlot3D. The mesh is obviously incorrect and some faces are turned inside out

region = RegionDifference[Cuboid[], Ball[{.5, .5, .5}, .2]];
RegionPlot3D[region, PlotStyle -> Directive[Opacity[0.8], FaceForm[Red, Blue]],
Mesh -> All, Lighting -> "Neutral", PlotPoints -> 13, MaxRecursion -> 0]


You can use BoundaryDiscretizeRegion instead

dr = BoundaryDiscretizeRegion[region];
Graphics3D[{EdgeForm[Opacity[0.1]], Opacity[0.5, Orange],
GraphicsComplex[MeshCoordinates[dr], MeshCells[dr, 2]]}]


• This does a good job when the size of the cavity is of larger, but when the cavity spheres are smaller, the cavity can sometimes not show up - or if it does it doesn't really look spherical. Is this a problem with the discretization process or with RegionDifference? – Jason B. Dec 14 '15 at 14:14
• @JasonB You can specify small MaxCellMeasure for BoundaryDiscretizeRegion. Unfortunately AccuracyGoal and PrecisionGoaldo not work here. – ybeltukov Dec 14 '15 at 14:24
• Oh~How nice job!Thanks very much.If you can extinguish the mesh of cell,I will appreciate you and accept it. – yode Dec 15 '15 at 1:50
• @yode - to get rid of the mesh, inside the Graphics3D, change it to say EdgeForm[Opacity[0.0]] – Jason B. Dec 15 '15 at 7:47
• Not sure why, when I change the EdgeForm[Opacity[...]] from 0.1 to 0.0, the mesh goes away. See the edits to my answer – Jason B. Dec 15 '15 at 9:09

Are you creating a region to do other calculations in, or are you just trying to plot a cube with several spherical cavities in it? If you are just trying to make the plot, I suggest not using RegionDifference. As I show below the plots, it is only an approximate function, and it fails a bit when the cavities are much smaller than the cube.

Here is the best plot I can make:

(* Define the cube, the cavities, and the difference region *)
length = 10;
o = {0, 0, 0};
dp = 0.25 (*m*);
rho = 2345 (*kg/m^3*);
rho = dp^3*rho;
cube = Cuboid[o, o + length];
cavity = Table[
Ball[RandomPoint[cube], RandomReal[{0.6 dp, 1.5 dp}]], 4];
region = Fold[RegionDifference, cube, cavity];

(* Then observe the difference between the plot with and without RegionDifference *)

Show[RegionPlot3D[#, PlotStyle -> Opacity[0.5], PlotPoints -> 20] & /@
Join[{cube}, cavity]]


versus

RegionPlot3D[region, PlotStyle -> Opacity[0.5], PlotPoints -> 20]


For comparison, let's look at the method given by @ybeltukov, with the mesh taken out as you requested,

dr = BoundaryDiscretizeRegion[region];
Graphics3D[{EdgeForm[Opacity[0.0]], Opacity[0.5, Orange],
GraphicsComplex[MeshCoordinates[dr], MeshCells[dr, 2]]}]


This is indeed better than just doing a RegionPlot3D on region, but it fails to make the cavities look spherical and it doesn't show all of the cavities. So in my opinion, the only way to make this plot well is to skip RegionDifference and just plot the cube and the spheres together.

But perhaps you are creating the RegionDifference in order to use it in calculations. One way to do this, is to do a three dimensional integral of some function over the region. Below, I run a test by integrating a simple function over the region defined by RegionDifference, and comparing this to the result you get by simply taking the difference of the integrals (integrate over the cube, then subtract the integrals over the cavities). The results should be identical:

NIntegrate[Sin[y + z - x], {x, y, z} ∈ region]
NIntegrate[Sin[y + z - x], {x, y, z} ∈ cube] -
Total@(NIntegrate[Sin[y + z - x], {x, y, z} ∈ #] & /@ cavity)
(* 6.57153 - 1.70511*10^-21 I *)
(* 6.5671 *)


They agree remarkably well, but not exactly. The integral over the RegionDifference took over ten times as long, and also returned convergence errors.

• Thanks for your pay attention to this problem.But as your see,the RegionPlot[Disk[]] should not give a right plot.BTW.what's meaning of your x^2 + Sqrt[y + z]? – yode Dec 15 '15 at 4:23
• @yode, I'm not sure what you mean by "should not give a right plot". Check the post, I updated it. I think the only way to get the plot you want is to not use RegionDifference – Jason B. Dec 15 '15 at 8:21
• I think your answer is nice,too.I'm feel so sorry.Thanks all the same. – yode Dec 15 '15 at 11:22
• @yode - nothing to be sorry about, whatever works for you – Jason B. Dec 15 '15 at 11:24
length = 10;
o = {0, 0, 0};
dp = 0.25 ;

region = Cuboid[o, o + length];
cavity = Table[Ball[{1, 1, 1}, RandomReal[{0.6 dp, 1.5 dp}]], {k, 4}];

Show[RegionPlot3D[region, PlotStyle -> Opacity[0.5], PlotPoints -> 10],
RegionPlot3D[cavity, PlotPoints -> 50, PlotStyle -> Opacity[0.7]]
] // Timing


You can change {1, 1, 1} with your definition of RandomPoint[].

• Thanks a lot,but the last region is expected,then I wanna visualize it.You just do the latter. – yode Dec 15 '15 at 0:07