# Inaccurate result from ImplicitRegion

I'm trying to describe the intersection of two cylinders by an ImplicitRegion (in order to use it in NDSolve later).

ImplicitRegion[
-0.5 <= x <= 0.5 && -0.5 <= y <= 0.5 && -0.5 <= z <= 0.5 &&
(x^2 + z^2 <= 0.04 || y^2 + (0.3 + z)^2 <= 0.04),
{x, y,z}];

Region[%]


The result is missing the sharp edges. How can I improve the accuracy of the intersection?

Please note that I get a simliar result when doing this using RegionUnion with Cylinder.

cyl1 = Cylinder[{{-0.5, 0, 0}, {0.5, 0, 0}}, 0.2];
cyl2 = Cylinder[{{0, -0.5, -0.3}, {0, 0.5, -0.3}}, 0.2];
RegionUnion[cyl1, cyl2];
ToBoundaryMesh[%];
MeshRegion[%]


Using discretized regions works fine for two cylinders, but when I add a 3rd cylinder I receive the error message:

The boundary surface is not closed because the edges Line[{{2120, 2126}, {2018, 2032}, {1919, 2018}, {1622, 2120}}] only come from a single face

cyl1 = Cylinder[{{-0.5, 0, 0}, {0.5, 0, 0}}, 0.2];
cyl2 = Cylinder[{{0, -0.5, -0.3}, {0, 0.5, -0.3}}, 0.2];
cyl3 = Cylinder[{{0, -0.5, 0.3}, {0, 0.5, 0.3}}, 0.2];
BooleanRegion[#1 || #2 || #3 &,
{BoundaryDiscretizeRegion[cyl1],
BoundaryDiscretizeRegion[cyl2],
BoundaryDiscretizeRegion[cyl3]}]


EDIT: Here is another example: If we just look at the surface mesh, we see, that the cylindrical faces are not connected, i.e. they do not form a connected region (even if I decrease the element size). So I cannot do an integral over that region.

ImplicitRegion[((x^2 + z^2 == 0.2^2)
|| ((y)^2 + (z + 0.3)^2 == 0.2^2)
|| ((y)^2 + (z - 0.3)^2 == 0.2^2))
&& ! (x^2 + z^2 < 0.2^2)
&& ! ((y)^2 + (z + 0.3)^2 < 0.2^2)
&& ! ((y)^2 + (z - 0.3)^2 < 0.2^2), {{x, -0.5, 0.5}, {y, -0.5,
0.5}, {z, -0.3, 0.3}}];
fluidSolidInterface =
ToBoundaryMesh[%, MaxCellMeasure -> {"Length" -> .1}]
%["Wireframe"]

• You could try easy one first and hard one later: BooleanRegion[#1 || #2 || #3 &, {BoundaryDiscretizeRegion[cyl2], BoundaryDiscretizeRegion[cyl3], BoundaryDiscretizeRegion[cyl1]}] – halmir Dec 27 '19 at 18:15

You can force Mathematica to use a finer mesh by using BoundaryDiscretizeRegion and giving the option MaxCellMeasure. Like so:

r =
ImplicitRegion[
-0.5 <= x <= 0.5 && -0.5 <= y <= 0.5 && -0.5 <= z <=
0.5 &&
(x^2 + z^2 <= 0.04 || y^2 + (0.3 + z)^2 <= 0.04),
{x, y, z}];
BoundaryDiscretizeRegion[r, MaxCellMeasure -> {"Length" -> .01}]


The smaller the value given for the length measure, the more accurate the result, but finer meshes come at the price of longer evaluation time.

I have created a GmshLink package as a workaround exactly for such questions. Please also see this answer for another nice example of use.

First we load the package and show path to directory containing GMSH executable.

Get["GmshLink"]
\$GmshDirectory = "path_to_directory\\gmsh-4.5.0-Windows64";


Define symbolic region and calculate its bounds (optionally).

reg = RegionUnion[
Cylinder[{{-0.5, 0, 0}, {0.5, 0, 0}}, 0.2],
Cylinder[{{0, -0.5, -0.3}, {0, 0.5, -0.3}}, 0.2]
]

bounds = RegionBounds[reg]


Create ElementMesh object with GmshGenerator function. It can accept different Options to further adjust the resulting mesh.

mesh = ToElementMesh[
reg,
bounds,
"BoundaryMeshGenerator" -> None,
"ElementMeshGenerator" -> {GmshGenerator},
MaxCellMeasure -> 0.05
]
(* ElementMesh[{{-0.5, 0.5}, {-0.5, 0.5}, {-0.5, 0.2}}, {TetrahedronElement["<" 9320">"]}] *)


We get a mesh of two cylinders with nice resolution of region intersection.

mesh["Wireframe"[
"MeshElement" -> "MeshElements",
"MeshElementStyle" -> FaceForm@LightBlue
]]


• Awesome! This is exactly what I'm looking for. Now I'm struggling with the ImportMesh package, which tells me that the GMSH format 4 is not supported. Any idea about that? – Mathias Luxner Dec 28 '19 at 7:59
• I'm glad that you found it useful. Currently ImportMesh can read only GMSH format 2 mesh files and GmshLink creates mesh files in that format. If you would like to import other GMSH mesh files, you should convert them to appropriate format first. Command line option for this is -format msh2. – Pinti Dec 30 '19 at 9:22

Boolean operations are easier to perform with discretized regions. Try this:

RegionUnion[
BoundaryDiscretizeRegion[cyl1],
BoundaryDiscretizeRegion[cyl2]
]
`