# How to make the boundary of a 3D region smooth?

I want to draw this region,but the surface is rough.I tried to find options to improve the surface but failed.

Region[RegionIntersection[Cylinder[{{0, 0, -2}, {0, 0, 2}}],
Cylinder[{{0, -2, 0}, {0, 2, 0}}],
Cylinder[{{-2, 0, 0}, {2, 0, 0}}]]]

which gives

besides, PlotPoints -> 200 seems have no place to put.

Sometimes, it is easier to discretize before intersecting.

BoundaryDiscretizeRegion@RegionIntersection[
Map[
BoundaryDiscretizeRegion[#,
MaxCellMeasure -> (1 -> 0.05)] &,
{Cylinder[{{0, 0, -2}, {0, 0, 2}}], Cylinder[{{0, -2, 0}, {0, 2, 0}}],
Cylinder[{{-2, 0, 0}, {2, 0, 0}}]}
]
]

Needs["NDSolveFEM"]
rr = RegionIntersection[Cylinder[{{0, 0, -2}, {0, 0, 2}}],
Cylinder[{{0, -2, 0}, {0, 2, 0}}],
Cylinder[{{-2, 0, 0}, {2, 0, 0}}]];
bmesh = ToBoundaryMesh[rr,
"ShapeSurfaceMeshOptions" -> {"LinearDeflection" -> 0.00125}}];
bmesh["Wireframe"[
"MeshElementStyle" -> Directive[FaceForm[Green], EdgeForm[]]]]

• Oh, that's great to know! I am eager to get my hand on the new version to try that! Commented Mar 20, 2020 at 19:10
• @HenrikSchumacher, it's bit like with TetGenLink and TriangleLink, OCCL forms a low level interface and with time more high level functions will make use of that. OpenCascade is a huge library and this first version scratches the surface, I think. It will take time to exploit more of it.... It's hopefully going to be easily accessible for you folks to make contributions too. Commented Mar 20, 2020 at 19:34
• Please comment on which OpenCascade app is needed to make this work. Commented Mar 20, 2020 at 21:56
• @murray, you just need Version 12.1 nothing else. ln the Documentation search for OpenCascadeLink for more information on this. Commented Mar 21, 2020 at 5:15
• @murray, it's a link to OpenCascade, does that help? The opencascade libraries now ship with mathematica 12.1 and the link, links to those shipped libraries. Commented Mar 23, 2020 at 5:20
ri = RegionIntersection[Cylinder[{{0, 0, -2}, {0, 0, 2}}],
Cylinder[{{0, -2, 0}, {0, 2, 0}}], Cylinder[{{-2, 0, 0}, {2, 0, 0}}]];

You can use PlotPoints as a suboption for Method options:

DiscretizeRegion[ri,
Method -> {"DualMarchingCubes", PlotPoints -> 150}]

In order to completeness,here we provide another ways to do this. Thanks @xzczd suggestion。

• MeshFunctions
Clear[f, g, h, fgh, ghf, hfg];
SetOptions[ContourPlot3D, PlotPoints -> 80, MaxRecursion -> 4,
Mesh -> {{0}, {0}, {0}}, BoundaryStyle -> None, Boxed -> False,
Axes -> False];
f = {x, y, z} |-> x^2 + y^2 - 1;
g = {x, y, z} |-> y^2 + z^2 - 1;
h = {x, y, z} |-> z^2 + x^2 - 1;
fgh = ContourPlot3D[
f[x, y, z] == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
MeshFunctions -> {g, h},
MeshShading -> {{Red, None}, {None, None}}, MeshStyle -> None];
ghf = ContourPlot3D[
g[x, y, z] == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
MeshFunctions -> {h, f},
MeshShading -> {{Yellow, None}, {None, None}}, MeshStyle -> None];
hfg = ContourPlot3D[
h[x, y, z] == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
MeshFunctions -> {f, g},
MeshShading -> {{Cyan, None}, {None, None}}, MeshStyle -> None];
Show[fgh, ghf, hfg]

• CSGRegion
CSGRegion["Intersection", {Cylinder[{{0, 0, -2}, {0, 0, 2}}],
Cylinder[{{0, -2, 0}, {0, 2, 0}}],
Cylinder[{{-2, 0, 0}, {2, 0, 0}}]}, BaseStyle -> Darker@Cyan]

• OpenCascadeLink (Thanks @user21)
RegionIntersection[Cylinder[{{0, 0, -1}, {0, 0, 1}}],
Cylinder[{{0, -1, 0}, {0, 1, 0}}],
Cylinder[{{-1, 0, 0}, {1, 0, 0}}]]];
"ShapeSurfaceMeshOptions" -> {"AngularDeflection" -> .01}] //
BoundaryMeshRegion;
RegionPlot3D[bm, ColorFunction -> "Rainbow", Boxed -> False]

Or

Needs["OpenCascadeLink"];
reg1 = Cylinder[{{0, 0, -2}, {0, 0, 2}}];
reg2 = Cylinder[{{0, -2, 0}, {0, 2, 0}}];
reg3 = Cylinder[{{-2, 0, 0}, {2, 0, 0}}];
"ShapeSurfaceMeshOptions" -> {"AngularDeflection" -> .01}] //
BoundaryMeshRegion;
bm // Volume;
RegionPlot3D[bm, ColorFunction -> "Rainbow", Boxed -> False]

Yet another way to discretize the Steinmetz solid is to directly derive the required inequalities using RegionMember[]:

ineq = Simplify[RegionMember[RegionIntersection[Cylinder[{{0, 0, -2}, {0, 0, 2}}],
Cylinder[{{0, -2, 0}, {0, 2, 0}}],
Cylinder[{{-2, 0, 0}, {2, 0, 0}}]],
{x, y, z}], {x, y, z} ∈ Reals]
-2 <= z <= 2 && x^2 + y^2 <= 1 && -2 <= y <= 2 && x^2 + z^2 <= 1 &&
-2 <= x <= 2 && y^2 + z^2 <= 1

which can then be fed to ImplicitRegion[]:

reg = ImplicitRegion[ineq, {x, y, z}];
BoundaryDiscretizeRegion[reg, MaxCellMeasure -> {"Length" -> 0.02}]

reg = Polygon[
Join[#, Reverse[#.DiagonalMatrix[{1, -1}]]] &@
Table[{x, Min[AngleVector[Mod[x, Pi/2]]]}, {x, 0., 2 Pi, 2 Pi/(25*8)}]]

ParametricPlot3D[{{Cos[t], Sin[t], h}, {h, Cos[t], Sin[t]}, {Cos[t], h, Sin[t]}},
Element[{t, h}, reg], ImageSize -> 600] // AbsoluteTiming

For version 13.x, you need to use VisualizationCore`ParametricPlot3D instead of ParametricPlot3D.